arXiv:hep-lat/0503040vl 30 Mar 2005 


HU-EP-05/14 

SFB/CPP-05-09 


Precision Lattice Computations in the 
Heavy Quark Sector 


DISSERTATION 

zur Erlangung des akademischen Grades 
doctor rerum naturalium 
(Dr. rer. nat.) 
im Fach Physik 

eingereicht an der 

Mathematisch-Naturwissenschaftlichen Fakultat I 
der Humboldt-Universitat zu Berlin 


von 

Herr Dipl.-Phys. Andreas Jiittner 
geboren am 19.08.1976 in Niirnberg 


Prasident der Humboldt-Universitat zu Berlin: 

Prof. Dr. Jurgen Mlynek 

Dekan der Mathematisch-Naturwissenschaftlichen Fakultat I: 
Prof. Thomas Buckhout, Ph.D. 

Gutachter: 

1. Dr. Jonathan Flynn 

2. Dr. Rainer Sommer 

3. Prof. Dr. U. Wolff 

eingereicht am: 26. Juli 2004 

Tag der miindlichen Priifung: 1. Oktober 2004 


Abstract 


The phenomenology of the pseudo scalar mesons Dg and Bg and of the vector 
mesons D* and B*, each of which contain a heavy and a light quark, was investi¬ 
gated in simulations of quenched lattice QCD. The work is particularly focused 
on the minimisation of all systematic errors within this approximation. 

The decay constants Fd^ and Fd* and the difference in the masses between the 
pseudo scalar Dg-meson and the corresponding vector meson D* were determined 
from the direct computer simulation of lattice QCD in large physical volume 
(L 1.5 fm). As an aside, the renormalisation group invariant charm quark 
mass Me could be obtained from the simulation results. 

A platform independent software was developed for the Monte-Carlo simula¬ 
tions of lattice QCD within the Schrodinger Functional. A number of simulations 
at different lattice constants allowed the extrapolation of the results to the con¬ 
tinuum. 

Since comparable simulations for the Bg- and the B*-meson are not feasible 
due to the large mass of the 5-meson, an interpolation in the meson mass to its 
physical point was carried out for the decay constant and the mass splitting. The 
interpolation was carried out between the static limit and the range of meson 
masses of order The desired observables were therefore determined and 

extrapolated to the continuum for altogether six meson masses. The functional 
form of the subsequent interpolation in the meson mass to the static limit was 
guided by the prediction of the Heavy Quark Effective Theory (HQET). In order 
to apply it to the results obtained in QCD, a set of conversion functions between 
HQET and QCD were derived and evaluated numerically with input from results 
in perturbation theory. 

The final results are Fd,, = 226(7)MeV, Fd* = 239(18)MeV, 
Fbs = 198(9)MeV, rriD* — tud, = 136(9)MeV, mej — = 63(6)MeV and 

Me = 1.60(3)GeV. The result for the renormalisation group invariant charm 
quark mass is equivalent to = 1.27(3)GeV. 

The analysis of the interpolation furthermore allowed to estimate, that the 
lowest order corrections to the static limit in HQET are relatively small. One 
therefore can expect HQET to offer a good approximation in the range of B- 
physics. 
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Zusammenfassung 


Die Phanomenologie der pseudoskalaren Mesonen Dg und Bg sowie der Vek- 
tormesonen D* und B*, welche jeweils ein schweres und ein leichtes Quark enthal- 
ten, wurde in numerischen Simulationen von Gitter-QCD unter Vernachlassigung 
virtueller Fermionschleifen untersucht. Besonderer Wert wurde auf die Kontrolle 
und Minimierung aller systematischen Fehler innerhalb dieser Naherung gelegt. 

Die Zerfallskonstanten Fd^ und F^j und die Massendifferenz zwischen dem Dg- 
und dem D*-Meson wurden aus der direkten Computersimulation von Gitter- 
QCD in groBem physikalischen Volumen (L 1.5 fm) bestimmt. Als Neben- 
produkt konnte auch ein praziser Wert der renormierungsgruppen-invarianten 
Charm-Quarkmasse ermittelt werden. 

Fiir die Monte-Carlo Simulationen von QCD auf dem Gitter, speziell im liier 
verwendeten Schrodinger Funktional, wurde eine plattformunabliangige Software 
entwickelt. Eine Reihe von Simulationen bei verschiedenen Gitterabstanden er- 
laubte die Extrapolation der Ergebnisse zum Kontinuum. 

Da vergleichbare Simulationen fiir das Bg- und B*-Meson aufgrund der groBen 
Masse des enthaltenen 6-Quarks nicht mbglich sind, wurde eine Interpolation in 
der Mesonmasse zu ihrem experimentell bekannten Punkt fiir die Zerfallskon- 
stante und fiir den Wert der Massendifferenz durchgefiihrt. Interpoliert wurde 
dazu zwischen dem statischen Limes (unendliche Mesonmasse) und dem Bere- 
ich von Mesonmassen in der GroBenordnung von mD^. Fiir insgesamt sechs 
Mesonmassen in diesem Bereich wurden die gewiinschten Observablen deshalb aus 
Simulationen von Gitter-QCD in groBem Volumen bestimmt und die Ergebnisse 
zum Kontinuum extrapoliert. Die Form der anschlieBenden Interpolation in der 
Mesonmasse zum statischen Limes wurde den Vorhersagen der Heavy Quark Ef¬ 
fective Theory (HQET) entsprechend gewahlt. Um diese auf QCD zu iibertragen, 
wurden Konversionsfunktionen zwischen HQET und QCD hergeleitet und mit 
Hilfe von Ergebnissen aus der Storungstheorie numerisch bestimmt. 

Die Endergebnisse sind FD^=226(7)MeV, FDj=239(18)MeV, Fb^= 197(9)MeV, 
mD* — ruDs = 136(9)MeV, me* — rriBs = 63(7)MeV und Me = 1.60(3)GeV. Das 
Ergebnis fiir die Quarkmasse ist Equivalent zu = 1.27(3)GeV. 

Aus der Analyse der so bestimmten Interpolationen lieB sich auBerdem ab- 
schatzen, daB die fiihrenden Korrekturen zum statischen Limes in der HQET 
relativ klein sind. Man erwartet deshalb, daB HQET im Bereich der R-Physik 
eine gute Naherung darstellt. 
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Chapter 1 
Introduction 


With the goal of a unihcation of quantum mechanics and special relativity, P. 
A. M. Dirac formulated his famous equation of motion for free spin-^ particles 
in 1928 [?]. The Dirac equation has solutions with both positive and negative 
energy. The latter found an interpretation in terms of anti-particles four years 
later, when C. D. Anderson discovered the positron in cosmic rays [?]. 

Today, based on Dirac’s hndings, elementary particles are described as the 
quanta of helds in local quantum held theories. The CPT-theorem [?,?,?], ac¬ 
cording to which the held theory has to be invariant under the combined ap¬ 
plication of charge conjugation (C), parity (P) and time-reversal (T), postulates 
an anti-particle to be associated to each particle. Leaving aside gravity due to 
its weak coupling to elementary particles at the energy scales accessible to ex¬ 
periments today, the electro-magnetic (QED), the weak and the strong (QCD) 
interactions have been combined in the Standard Model of elementary particles 
which has an underlying local SU(3)c x SU(2)l x U(1)y gauge-symmetry. The 
gluons, the gauge bosons in QCD, and the photons, the gauge bosons in QED, 
couple to both the left-handed and the right-handed fermions. C and P are 
therefore good quantum numbers in these cases. In contrast, the gauge bosons 
and Z of the electro-weak sector only couple to the left handed fermions and 
therefore parity is violated. But at least the combination CP for the time being 
seemed to be a symmetry of the electro-weak interactions. 

However, the Standard Model for three generations of quarks 



comprises the possibility of CP-violation. Although not conhrmed experimen¬ 
tally, the generation of particle masses in the Standard Model is explained by 
the Higgs mechanism. The scalar Higgs held interacts with the lepton and the 
quark helds through Yukawa couplings. By spontaneous electro-weak symmetry 
breaking, the Higgs held acquires a non-vanishing vacuum expectation value and 
thereby dynamically generates a mass term for all helds. The physical quark 
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fields M, d, s, c, 6, are in the mass eigenstate basis, where the associated mass 
matrix is diagonal. Their relation to the quark states in the weak eigenstate ba¬ 
sis is given by a 3 x 3 unitary matrix, the Cabbibo-Kobayashi-Maskawa (CKM) 
mixing matrix [?] 



to the freedom to redehne the phases of the quark mass eigenstates, leaving a 
single physical phase (5 km, the Kobayashi-Maskawa phase. In the case of the CP- 
invariance of the electro-weak interactions, the phase has to vanish. However, 
in 1964, Cronin, Fitch and Christenson [?] found experimental evidence for a 
CF-violating 27r-decay of the neutral Kaon. 

Since the discovery of CF-violation, a lot of effort has been put into preci¬ 
sion measurements of the CKM-matrix. New experiments, like CLEO-c [?] and 
the B-factories BaBar [?] and BELLE [?] have been set up for high precision 
measurements of CF-violating effects. 

The question arises as to why these measurements are interesting. 

Once the parameters of the Standard Model have been hxed by experiment, 
the consistency of the theory can be checked. In particular, as CF-violation in the 
Standard Model is induced only by the phase (5 km, hs measurement in one process 
will constrain the CF-violation allowed in other processes. For example, the 
CF-violation in the decay B —»• 'ijjKs is related to the CF-violation in K —>■ ttz/P 
[?]. If inconsistencies were discovered by experiment, this would be a sign for 
physics beyond the Standard Model. Indeed, super-symmetric models predict 
CF-violating effects, that would exceed the magnitude of CF-violation allowed 
by the Standard Model [?]. 

According to one of Sakharov’s three conditions [?], all of which must be met 
in order to allow for baryogenesis, CF must be violated in order to favor baryon 
over anti-baryon production. Therefore, a deeper understanding of CF-violation 
in connection with electro-weak symmetry breaking in the early universe may 
help in understanding the observed asymmetry between matter and anti-matter 
in the observable part of the universe. 

In order to assess possible inconsistencies in the flavor physics of the Standard 
Model, precision measurements on the one hand and precision predictions from 
theory on the other hand are required to reduce the error on the elements of the 
CKM-matrix. For example, the CLEO-c experiment intends to measure precisely 
the branching ratio of the leptonic decay of a Dg-meson, 



(1.3) 


In this expression, which is correct up to radiative corrections, Cp is the Fermi 
constant, Ed,, td, and m-D^ the lifetime, the decay constant and the mass of 
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the Ds-meson respectively and mi is the mass of the lepton taking part in the 
decay. In order to extract the valne of the matrix element |14s|, theorists have 
to deliver precise predictions for the decay constant Fd^. However, the Dg-meson 
is particular, as iKsI can also be obtained very precisely from constraints on the 
CKM-matrix [?] and therefore the experiment can also measure Fd^. Hence, its 
study offers the possibility to test precision predictions from the lattice and to 
ensure, that the same techniques applied to experimentally less explored meson 
systems produce reliable results. 

In similar ways, from the measurements of the mixing Bd <->■ Bd and Bg ^ Bg, 
the product of CKM-matrix elements | ltd 11 Kb I and the ratio |Kd|/|Ks| could be 
determined. To do this, theorists would have to predict accurately the products 
of the decay constant and the square root of the bag-parameter, FB^A/Bd and 
Fbsa/B^ from the study of semi-leptonic decays [?]. 

On the theory side, in order to determine these low-energy quantities, QCD 
sum rules [?,?] or relativistic quark models [?] may be applied. Lattice QCD 
however, where a Euclidean space-time lattice is used as the regulator for QCD, 
appears to be the most promising approach. Physical quantities can be obtained 
as the expectation values of observables evaluated on an ensemble of held con- 
hgurations, which have to be computer-generated by means of a Monte-Carlo 
simulation [?,?]. 

For precision lattice phenomenology, a number of systematic uncertainties 
present in the lattice approach have to be taken into account [?]: 

• Finite volume effects - Simulations of lattice QCD are numerically very 
costly and therefore the size of the lattices which can be simulated is limited. 
Especially for light quarks, where the Compton wavelength is large, one has 
to make sure that results are not affected by the presence of a space-time 
boundary. 

• Continuum limit and cutoff effects - Lattice QCD has to be simulated at 
various hnite lattice spacings in order to allow for a controlled extrapolation 
to the continuum limit. This has to be done along a line of constant physics, 
which amounts to renormalizing the theory. Systematic effects can be re¬ 
duced by renormalizing the theory non-perturbatively [?]. Furthermore, 
the approach to the continuum can be accelerated and systematic effects 
can be reduced by non-perturbatively improving the theory [?]. In this way, 
lattice artifacts vanish quadratically in the lattice spacing [?, ?] instead of 
linearly, as it is the case for example for standard Wilson fermions [?]. 

• Extrapolations to physical quark masses - Lattice computations become very 
costly for light quarks {u and d) and for heavy quarks with masses above the 
charm quark mass. Light quarks require a large physical volume and at the 
same time, they are particularly costly with current algorithms due to the 
required inversions of the badly conditioned Dirac operator. For the heavy 
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quarks, mass dependent cutoff effects at finite lattice spacing are a source 
of concern due to their short Compton wave length. Both problems can 
be circumvented by simulating at unphysical but less costly quark masses 
and then extrapolating the results to the physical point. The form of the 
extrapolation is suggested by effective theories. Chiral perturbation theory 
is commonly used for extrapolations to light quark masses, whereas heavy 
quark effective theory (HQET) suggests a polynomial expansion of observ¬ 
ables in the inverse quark mass for heavy quarks. The extrapolations in 
each case have to be done carefully to avoid uncontrolled systematic errors. 

• Excited states - In lattice computations, physical observables are mostly 
extracted from the time dependence of correlation functions. One does 
not know exactly how to construct particle wave functions in QCD which 
would allow for a projection onto particular states of the spectrum. There¬ 
fore, contributions of excited states to the desired observables cannot be 
completely avoided. The magnitude of such contributions can however in 
some cases be estimated and controlled by an accurate data analysis. 

• Quenching - The generation of a representative ensemble of gauge conhgu- 
rations in the lattice Monte-Carlo simulation is very costly in the case of full 
QCD. Full QCD means that the simulation takes into account both virtual 
gluon loops and virtual quark loops. In the quenched approximation, one 
neglects the contributions of virtual quark loops. This reduces the costs 
immensely, at the expense of signihcance of the results: Quenched QCD is 
an uncontrolled approximation to QCD and precise and reliable results for 
phenomenology cannot be obtained. But still, simulations in the quenched 
theory give estimates for phenomenological quantities which in some cases 
are surprisingly good and are an important tool to assess techniques for 
later use in the full theory. 

In this work, a feasibility study of precision lattice computations in the heavy- 
fiavor sector of quenched QCD has been accomplished. A particular emphasis 
was placed on keeping the above sources of systematic errors, apart from quench¬ 
ing, under control. The focus of the study was on the phenomenology of the 
heavy-light pseudo scalar mesons Dg and Bg and the vector mesons D* and B*. 
In particular, the simulations aimed at the computation of the leptonic decay 
constants Fd^, Ed*, Fb^ and Fb*, which are important for the determination of 
CKM-matrix elements. 

Simulations of lattice QCD for heavy-light mesons in large physical volume 
(L > 1.5 fm) and with a controlled continuum extrapolation currently can be 
accomplished at the physical point of charm and strange as the heavy and light 
quark respectively. Thus, the Dg-meson can be simulated directly and systematic 
effects due to extrapolations in the quark mass can be avoided. Since very precise 
data for the matrix element iKsI exists, CLFO-c will determine Fd^ to a precision 
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Figure 1.1: Using data from lattice simulations of QCD with heavy quark masses 
around the charm quark mass in the continuum limit (bold dashed line with 
error band), together with predictions or simulation results in the static limit 
rriQ —>■ oo, allows for an interpolation to the physical point of the 6-quark mb- 


below 2% until 2005 and thereby offers a very accurate test of the lattice approach 
and all the applied techniques. 

Once one has gained conhdence in the lattice computations, the techniques 
can be applied to sectors of the Standard Model, like the Bg-mesons, which are 
not easily accessed through experiments. The large mass of the 6-quark however 
does not allow for a direct lattice simulation of the Bg mesons. Instead, the 
following procedure, which is sketched in figure 1.1, can be applied. In addition 
to the simulations at the physical point of the Dg-meson, one simulates for the 
desired meson observable also at a number of unphysical heavy quark masses 
around charm (indicated by the dashed bold line). After the continuum limit 
has been taken, the data would allow for an extrapolation to the physical point 
of the 6-quark mass. However, the range of heavy quark masses accessible to 
relativistic simulations of QCD is limited. Therefore, such an extrapolation has 
little signihcance and systematic effects cannot be estimated. 

Fortunately, the extrapolation can be constrained further. HQET makes ex¬ 
act predictions for some mesonic observables in the limit of inhnite heavy quark 
mass, the static limit. Also, results for the decay constant in the continuum, 
obtained from lattice simulations in the static approximation with reasonable 
statistical errors, exist [?]. Furthermore, HQET predicts the mass dependence of 
mesonic observables as a polynomial in the inverse heavy quark mass. Combining 
predictions from HQET in the static limit with lattice QCD simulations in the 
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charm region allows for a controlled interpolation to the beauty region [?,?,?], 
while keeping systematic effects under control. Furthermore, in assessing the 
functional form of the resulting interpolation, an estimate of the order of mag¬ 
nitude of the leading order coefficients in the heavy quark expansion is possible 
and at the same time constitutes a test for HQET. 

HQET [?,?,?] is an effective theory for QCD with heavy quarks of mass mq S> 
Aqcd- a short motivation and its definition are given in chapter 2. Conversion 
functions that relate observables of heavy-light mesons in HQET in the continuum 
to their analog in QCD have been derived, computed and parameterized in terms 
of the renormalization group invariant heavy quark mass Mq. These conversion 
functions will hnally enable an interpolation between results from QCD in the 
charm region and the static limit guided by the predictions from HQET. 

Section 3 establishes the QCD Schrodinger Functional as the preferred frame¬ 
work for the lattice computation of decay constants, meson masses and quark 
masses, for which explicit expressions in terms of correlation functions and £- 
nally in terms of quark propagators will be derived. 

A major part of this work consisted in obtaining a platform independent 
program code that can accomplish all the necessary computations. The program 
code presented in chapter 4 is based on the MILC collaboration’s lattice gauge 
theory code [?]. All major changes, improvements and tests of the code will be 
discussed. 

Chapter 5 discusses and tabulates all parameters for the Monte-Carlo simu¬ 
lations. 

The analysis of all data is described in chapter 6. After the discussion and 
estimation of all sources of systematic errors, the continuum extrapolation for 
the Dg- and the D*-meson are presented. Afterwards, the results from the inter¬ 
polation to the static limit are discussed. On the one hand, the values of meson 
observables at the physical point of the Bg have been obtained and on the other 
hand an estimate for the order of magnitude of the leading order coefficient of 
the heavy quark expansion will be given. All results are discussed and compared 
with other lattice computations and experiment. 

The last chapter summarizes all findings of this work and gives an outlook. 



Chapter 2 

Non-perturbative test of HQET 
with QCD 


After a short motivation of HQET in section one, section two describes the com¬ 
putation of conversion functions, which are necessary to relate matrix elements in 
QCD to those in HQET. They will allow for an interpolation in the heavy quark 
mass between the region of the charm quark mass and the static limit, guided by 
predictions from HQET. The necessary relations are given in section three. 


2.1 Heavy quark effective theory 

The typical energy carried by the light constituents in mesons (u-, d-, s-quarks 
and anti quarks and gluons) is of order Aqcd ~ 200 MeV. The phenomenology 
of mesons containing a light quark q and a heavy quark Q with tuq ^ Aqcd^ 
(cf. hgure 2.1) is therefore governed by the two different energy scales mg and 
Aqcd- With the heavy quark’s Compton wavelength being Xq ~ l/'m-g, the 
gluons cannot resolve the heavy quark’s quantum numbers - the light degrees of 
freedom are blind to spin and flavor (mass) of the heavy quark, leading to heavy 
quark spin and flavor symmetry. For example the experimentally determined 
spin splittings [?] 

m|. — m| 0.49 GeV^, . . 

~ 0.55 GeV^, 

and the mass splittings [?] 

rriB, - = (90 ± 3) MeV, . . 

m’Da — = (99 ± 1) MeV, 

for different heavy-light mesons are approximately the same. One expects these 

particular choice for the quark mass definition will be done in section 2.2.2. At this point 
tuq may for example be the heavy quark’s pole mass. 


7 



CHAPTER 2. NON-PERTURBATIVE TEST OF HQET WITH QCD 


u/d 



j I..I 


10 



J_I. . I _LJ_I_I. . I _ I _I_I_ I I I 

100 1000 m^(2 GeV)/MeV 


Figure 2.1: Quark mass ranges in the MS-scheme of dimensional regularization [?]. 

symmetries to be exact for heavy-light mesons with one inhnitely heavy, or static, 
quark. The symmetry breaking which can be observed experimentally at a hnite 
but large heavy quark mass can be interpreted as the consequence of small per¬ 
turbations to the theory with a static quark due to the interaction with the 
chromo-magnetic and chromo-electric helds mediated by soft gluons. This idea 
has been formulated in terms of an effective theory, the heavy quark effective 
theory (HQET) [?, ?, ?] which will be derived briefly in the following. For an 
extensive derivation the reader may refer to the reviews [?,?] and [?]. 

Restricting the study on heavy-light mesons with momentum p, containing 
one flavor of heavy quarks Q{x) and one flavor of light quarks q{x), the starting 
point is the QCD path integral 

is the SU(3) Yang-Mills Lagrangian, and is the QCD Lagrangian for 
quark helds coupled to the gauge held U in the adjoint representation, 

D^['0(a:),'0(a;), f/(a;)] = ip(x')(iI/> + m)ip(x'), (2.4) 

with the Dirac operator .0. The focus will now be on the heavy quark Lagrangian. 
The heavy quark in the meson is approximately on shell and therefore behaves like 
a free particle moving at four-velocity v. Removing the space time dependence 
of a solution of the free Dirac equation, the four-component Dirac held Q(x) can 
be rewritten in terms of the large and small component helds^ hy(x) and Hy{x) 
by 

hy{x) = e^^Q^-^PlQ{x) andHy{x) = P^Q^-^PlQ{x). (2.5) 

P^ and P^ are the projection operators 

1 zb 

_ n = (2.6) 

^This nomenclature stems from the free Dirac theory, where in the non-relativistic limit 
E —> mP, the upper components of the Dirac spinor remain of 0(1) while the lower components 
vanish. One therefore refers to the upper components as the “large components” and to the 
lower components as the “small components”. 
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Xo 


Figure 2.2: Virtual fluctuation of a heavy quark. 


The time dependence of the helds h{x) and H{x) is then expected to be deter¬ 
mined by the residual momentum k = p — mqv which is of order Aqcd- The 
heavy quark will only be considered in its rest frame throughout this work and 
therefore = (1,0, 0,0). In this case, h{x) = hy{x) corresponds to the upper 
components of Q{x) and H{x) = H^{x) to the lower components. 

The small components H{x) of the heavy quark held Q{x) only become rele¬ 
vant at high energies and are the origin of the short distance effects - for example, 
effects involving pair creation of heavy quarks or the zig-zag depicted in hgure 
2 . 2 , where the intermediate state has an energy that differs from the initial one 
by at least 2mQ and therefore propagates only over a short distance. 

In terms of the helds H{x) and h{x), the heavy quark Lagrangian in the rest 
frame can be rewritten as 


C[Q{x),Q{x),U{x)] = h{x)iDoh{x) — H{x){iDo + 2mQ)H{x) 


+ h{x)ilp±H{x) + H{x)ilj)±h{x). 


(2,7) 


with D± = {0,D). The large component helds h{x) do no longer have a mass 
term, whereas the small component helds H{x) appear with a mass term with 
twice the heavy quark mass. It is this term which will be eliminated in the 
construction of the ehective theory. 

By Gaussian integration, which in this case is equivalent to applying the 
classical equation of motion 


(iDo + 2mQ)H[x) 


( 2 . 8 ) 


the small component helds can be eliminated and one arrives at the non-local 
ehective Lagrangian 

C^^[h{x),h{x),U{x)] = h{x)iDoh{x)+ h{x)ilp ±-— ^ iuo ilp±h{x). (' 29 ) 

The second term in this Lagrangian represents the virtual processes suppressed by 
at least l/2mQ. In momentum space the operator that acts on h{x) corresponds 
to powers of the momentum. As the residual momenta of the heavy quark held 
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h{x) are small with respect to the heavy quark’s mass, the quotient in the second 
term can be expanded in Dq/ttiq by means of a derivative expansion which results 
in a an effective Lagrangian, in which the operators are ordered in powers of 
1/mg. This is the HQET Lagrangian. Up to the 1st order in l/mg it reads^ 

[h{x), h{x), U(x)] = C^''^^[h{x),h{x),U{x)] 

( 2 , 10 ) 

+ l‘{x),U(x)] + 0(l/ml) 

with 


C^^^^[h{x),h{x),U{x)]= Osti,t[Hx),h{x),U{x)] = h{x)iDoh{x), 
UjmQ [Kx), h{x),u{x)]=0\f^^ [h(a:), h{x),U{x)]+Off^^ [h(x), h{x),U{x)] 

= h{x)iD‘^h{x) + h{x)iS ■ B{x)h{x). 


( 2 . 11 ) 

The S'* are the generators of spin SU(2) rotations and can be chosen as 

= 1 ^ S^] = (2.12) 

where the (Tj are the Pauli matrices (cf. appendix A). B^{x) = — |e*-^^G*-^((r) is 
the chromo-magnetic gluon held where [iD°‘{x),iD^{x)] = igG^^{x) is the gluon 
held strength tensor. The term with is responsible for huctuations of order 

Aqcd in the heavy quark’s motion and describes the coupling of the heavy 

quark’s spin to the chromo-magnetic held. Both terms introduce the leading 
order havor and spin symmetry breaking interactions at hnite heavy quark mass, 
which were mentioned at the beginning of this chapter. 

The theory with the Lagrangian (2.10) is not renormalizable by a hnite num¬ 
ber of counter terms. Due to the presence of couplings with negative mass dimen¬ 
sion, terms of a given order in l/mg may mix with terms of higher order under 
renormalization [?] and an inhnite number of counter terms would be necessary. 

Thus, one expands the Boltzmann-factor in the corresponding path integral 
in the heavy quark mass l/mg, 

_ f p—i J d‘^xC’^^‘^^[h(x),h(x),U(x)] 

Jh,h 

^ S<'-''xCii„„[h(x),h(x),U{x)\ + 0(llw?g)^ . 

_ (2.13) 

•^Higher order terms will not be considered in this thesis. 
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From power counting one concludes, that the static theory dehned by is 
renormalizable with a finite number of parameters. 

In the same way as for the derivation of the HQET Lagrangian, an operator 
0^(x) containing heavy quark degrees of freedom, at tree-level can be expanded 
in a power series in l/mg, 

0^(x) = Of (x) + + 0(l/m|), (2.14) 

Q 


This may for example be done for the heavy-light axial vector current A^(x) = 
9(3^)7 m 75Q(^) (X = PS'^) and the vector current Vfj,(x) = q{x)'jfj,Q{x) (X = V) 
which then at leading order are defined as 


OPS(x) = A^{x) = q{x)-f^'j5h{x), 
0^{x) = V^{x) = q{x)-f^h{x). 


(2.15) 


Unlike the analog weak current operators in QCD, Oq^^{x) and become 

scale dependent under renormalization. Also the chromo-electric moment 
receives a scale dependence. In contrast, stays scale independent due 

to re-parameterization invariance [?,?]. 

For the cases X=PS, V and spin one then writes 


Ol{x,^^) = Z^{^^)0^{xl (2.16) 

with the renormalization constant whose scale dependence is determined 

by the renormalization group equation 

(2.17) 

d/i 

The renormalized coupling g{fi) is the one in the MS-scheme of dimensional reg¬ 
ularization and the anomalous dimension 7 ^’^®(( 7 ) has the generic perturbative 
expansion 

7 ^’^(d) = -7^^^ - 7^^" - 7pd® + • • • • (2.18) 

It is equivalent for X=PS and V and has been determined in the MS-scheme of 
dimensional regularization at one-loop in [?,?], at two-loop in [?,?] and at three- 
loop precision in [?]. For X = spin, the one-loop anomalous dimension is given 
in [?,?] and at two-loop in [?,?]. The corresponding coefficients are given in table 
2.1. The vacuum expectation value of an operator 0^(x,/u) in HQET then takes 

^This common notation refers to the transformation properties of A Ax) under parity (odd), 
which are the same as for a pseudo scalar. 






12 


CHAPTER 2. NON-PERTURBATIVE TEST OF HQET WITH QCD 



X=PS,V 

X=spin 

X,MS 

7o 

1 

6 

(47r^) 

(47r)^ 

X,MS 

7i 

(254 . 56 \ 1 

V 9 ^ 27n^) (47r)'‘ 

68 

X,MS 

72 

12.941 

{4^2)3 



Table 2.1: Coefficients for the 3- resp. 2-loop anomalous dimension for renor¬ 
malized heavy-light quark currents (axial vector and vector current) and the 
chromo-magnetic moment of a heavy quark. 


the form 

+ Jd^y f^)) ) 

+ 0{l/ml), 

(2.19) 

where the operator expectation values have to be understood in the theory dehned 
by the path integral 

Z = J- _ ^ ^ (2.20) 

2.2 Matching the effective theory to QCD 

By explicitly integrating out the short distance physics associated with the heavy 
quark in the last section, an effective theory for heavy quarks has been derived 
which one expects to correctly describe the long-distance physics of QCD. 

It is known from QCD, that quarks couple to gluons which can have virtual 
momenta as high as the quark mass. In HQET, when taking the limit mq —>■ oo, 
this introduces logarithmic divergences for example in weak matrix elements. 
Those matrix elements therefore have to be renormalized. 

The matching of the effective theory to QCD amounts to reintroduce the high 
energy behavior of matrix elements in HQET in terms of Wilson coefficients. 
They allow to dehne conversion functions C, which relate QCD matrix elements 
for heavy quarks of mass mq to the corresponding renormalization group invariant 
matrix elements in HQET. 
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2.2.1 The conversion functions CxirriQ) for 

X = PS, V, PS/V, spin 

The Wilson coefficients are defined by the relation between the matrix element 
of the corresponding operator Oji{x,mQ) in QCD, containing heavy degrees of 
freedom of mass mg, and the operator in the effective theory, renormalized at 
the scale p, 


(0^{x,mQ)\ = C'x(mQ,/i)/(0^)R(x,/i)\ 

\ / QCD \ / stat 

/ / stat I 

( 2 . 21 ) 

The coefficient Bx{mQ,fi) is mentioned for completeness but will be of no rele¬ 
vance for this work. In practice, one determines the Wilson coefficients C'x(mQ, p) 
in perturbation theory at the scale p = mg from a comparison or matching of 
suitable matrix elements in the full and in the effective theory^. Here, mg is the 
heavy quark’s pole mass. As the pole mass does not have a well dehned perturba¬ 
tive expansion [?], it will be replaced by the renormalization scheme independent 
renormalization group invariant quark mass Mg in the next section. 

CxijnQ^mo) depends on the particular Dirac structure of the operator 
Or(x, /i) and has been determined in perturbation theory for a number of heavy- 
light current matrix elements. The coefficients have an expansion in a power 
series in the renormalized coupling 

C'x(mg, mg) = 1 + c^f{mQ) + c^^^(mg) + ... . (2.22) 

For the axial vector and the vector current, the one-loop computation has been 
accomplished in [?] and at two-loop precision it is given in [?]. In the case of the 
kinetic term, Ckin(p-,pO = 1 holds due to re-parameterization invariance [?,?]. 
For Cspin(mg,mg) only the one-loop coefficient is known [?]. The factors cf and 
cf for the quenched theory are collected in table 2.2 for the phenomenologically 
important cases X=PS, V and spin. 

The scale dependence of the Wilson coefficients derives from the renormaliza¬ 
tion of the associated heavy quark current (cf. section 2.2). After integrating the 


^In [?], a method, how to do the matching non-perturbatively has been suggested. 
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X 

matrix element 


„x 

PS 

^'ps = (0|g7o75Q|PS) 

2 1 

3 477^ 

AO 1 

(4^2)2 

V 

$v = (0|g7oQ|V) 

4 1 

3 477^ 

115 ^ 

rr-o (4^2)2 

spin 

$"pP“= (PSIQ5-.BQIPS) 
(V|Q5-BQ|V) 

13 1 

6 477^ 

— 


Table 2.2: Coefficients for the matching factors Cx(mQ,mQ) in the qnenched 
theory. 

renormalization gronp eqnation (2.17) one gets the relation 

{ ff(0 

I 

Here, 7^’^®(s') is the anomalons dimension introdnced in the last section and P{g) 
is the anomalous dimension of the renormalized coupling g{fi) in the MS-scheme 
of dimensional regularization which is known at 4-loop accuracy [?], 

Pig) = -bo9^ - hg^ - b2g'^ - hg^ - . (2.24) 

The leading coefficients are bo = ll/(47r)^ and fei = 102/(47r)^ and the higher 
order coefficients are collected in appendix B. 

To eliminate any dependence on the renormalization scale in the relation 
between matrix elements in QCD and in HQET, it is convenient to take the limit 
yU —>■ oo in the above expressions. The Wilson coefficients then relate matrix 
elements in QCD to the renormalization group invariants 

0 ^ 01 ( 2 ;)= lim (2.25) 

//—>00 L J 

in HQET and one can write 

Cximq, n)0^{x, /i) ^ C'x(mQ)ORGi(a:). (2.26) 

2.2.2 Computation of Cx{Mq/Aqcb) 

Since the pole mass mg has a badly behaved perturbative expansion due to 
non-perturbative infrared effects [?], it will now be eliminated in favor of the 
dimensionless ratio between the renormalization group invariant quark mass Mg 
and the Agco-parameter as the new argument of the conversion functions®. Mg 

®Since only the case Nj = 0 was considered in this thesis, the non-perturbatively determined 
value Aqcd = = 238(19) MeV (quenched) [?] is used. 
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X,MS 


ig) 


Pig) 


(2.23) 
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is scale- and scheme-independent. It is defined via the limiting behavior of any 
renormalized mass 

Mq = hm {[ 26 o^ 2 ^/i)]-"°/( 2 feo)^(^)| ^ (2.27) 

/i—^■OO 

where do = 8/(47r)^ is the universal leading order coefficient of any quark mass 
anomalous dimension. How rriQ and Mq are related to each other in detail is 
explained in the appendix B. 

As an intermediate step, in the computation of the coefficients C'x(Mq/Aqcd), 
one defines the conversion functions parameterized with the renormalized mass 
m* = m(m*) in the MS-scheme, 






(2.28) 


The anomalous dimension 13{g) and the the anomalous dimension for X = PS, V, 


7^(f/) 


-7o"/ 






(2.29) 


will always be taken at 4- respectively 3-loop precision. The difference to taking 
the 3-loop /3-function instead, turned out to be tiny. The perturbative error 
introduced by 7^((?) was estimated with half the difference between the values 
for Cx obtained with the 2-loop and the 3-loop expression. For X=PS, V, the y* 
are dehned as 


7o 

X,MS 

= 7o , 


7f 

= 7f’^ + 26ocf, 

(2.30) 


= 72^’^ + 45o(4 + 7ok) + 2hcf - 2bo[cfr. 



All the coefficients are collected in the tables 2.1 and 2.2. 7 ^( 5 ') contains a 
contribution which has been derived from the matching (2.22) of the HQET 
operators and a contribution that originates from a re-parameterization: The 
matching was originally done at the matching scale given in terms of the heavy 
quark’s pole mass mq. Using the ratio mQ/Wi^, which is known at three-loop 
precision [?,?,?] (cf. appendix B), the pole mass can be replaced by m*. However, 
given the anomalous dimensions of the currents to three-loop order, only the one- 
loop term actually contributes to 7 ^( 5 ') (appearing as the piece proportional to 
k = —l/(37r^) in equation (2.30)). 

The chromo-magnetic operator Q{x)S ■ B{x)Q{x) in the heavy quark expan¬ 
sion is multiplied by the inverse pole mass. Since the preferred expansion param¬ 
eter for HQET in this thesis is Mq rather than mq, the corresponding conversion 
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function Cgpin(MQ/AQCD) must also include the factors m^/niQ and Mq/iti^: in or¬ 
der to cancel the factor l/mq in favor of 1/Mq. Using the relation (cf. appendix 

B) 


MQ/m^ = [2bog‘^{m,)] exp 




( 2 . 31 ) 


where _ 

= -g'^do - g'^di + ... (2.32) 

denotes the quark mass anomalous dimension in the MS scheme in QCD known 
up to four-loop precision [?,?], one then obtains for X=spin 


spin spin,MS i 

7o = 7o -«o, 


spin spin,MS 

7i = 7i 


di 2bo{cf^^ k), 


(2.33) 


where di = 404/(3(47r)^). For the case X = PS/V, all but the contributions from 
the matching cancel and one gets 


Cps/v(mQ = exp 


g[m*) 

I dp 


0 


7^^(^)-7^(l^) l 

(d{g) j ' 


(2.34) 


Using (2.31), one hnally changes the argument of the various C'x(m*) to the 
renormalization group invariant ratio Mq/Aqcd and arrives at expressions for 
the conversion functions 


Cx(Mq/Aqcd) with X = PS, V, PS/V and spin. (2.35) 


For practical purposes, such as repeated use in the hts of the heavy quark mass 
dependence of QCD observables that will be considered, a parameterization of 
all conversion functions in terms of the variable 

1 

In (Mq/Aqcd) 

was determined from a numerical evaluation. This parameterization is suggested 
by the asymptotic behavior of the conversion functions 


(2.36) 


Cx(M,/Apco) {l + O } 

(2.37) 

for X=PS, V and spin. The functions decompose into a prefactor encoding the 
leading asymptotics as x —0, multiplied by a polynomial of appropriate order 
in X. The results are given in table 2.3 and plotted in hgure 2.3. 
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Figure 2.3: Plots of the perturbative matching coefficients Cx, X=PS, V, PS/V 
and spin for 1-Ioop, 2-loop and 3-loop 7 -function (dotted, dashed and solid line 
respectively). 


The parameterizations of the matching factors deviate little from the numer¬ 
ical data and the error from perturbation theory is not too large (cf. table 2.3). 
Taking as an estimate for it half the difference between the numerical data based 
on the n-loop 7 -function and the (n — l)-loop 7 -function, the functions Cps, Cy, 
Cps/y have an error in the range of 1 % — 5.3%. 

The results for Cgpin with the 1-loop and the 2-loop anomalous dimension 
differ only little. As this may be accidental, the size of the three-loop term in 
Cps has been taken as the uncertainty in order to arrive at a conservative estimate 
for the error. A better error estimate would require the knowledge of 72 ^™. 

2.3 Combining QCD and HQET 

The decay constants for pseudo scalar mesons Fps(mps) and vector mesons 
Fv(mv) are dehned as 


(0|A^(0)|P) = ipf^Fps, 

(0|V(0)|V(A)) = ze/mvFv, 

where |P) and |V) are zero-momentum states with the quantum numbers of a 
pseudo scalar and a vector meson, respectively, and the | 0 ) represents the ground 
state of the gluonic vacuum. In the case of the vector meson decay constant, 
is a polarization vector, mps and my are the associated meson masses. With 
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parameterization of matching factors 

precision 

A max 
"^num. 

A max 
^pert. 

Am — 

Am — 

Cps(x) 

= x'^o®/(2^>o)(i _ 0.068X - 0.087x2 q.OTTOx^) 

3-loop 7 ^® 

10-2% 

1.8% 

0.7% 

1.1% 

C'ps(x) 

= xT'o®/(2feo)(l- 0.065x + 0.048x2) 

2-loop 7 ^® 





C'v(x) 

= x^o /(2bo)(i _ 0.196X - 0.222x2 + O.lOSx^) 

3-loop 7 ^ 

10-2% 

5.3% 

1.9% 

3.0% 

Cvix) 

= xT'^/(2fco)(i_o.i80x +0.099x2) 

2-loop 7 ^ 





Cps/y{x) 

= 1 + 0.124X +0.187x2 -0.102x3 

2-loop matching 

10-2% 

2.9% 

1.1% 

1.8% 

C'ps/v(^) 

= 1 + 0.117X - 0.043x2 

1-loop matching 





C^spin(^) 

= xT'r“/(2bo)('i + o.087x-0.021x2) 

2-loop 7 "P“ 

0.3% 

- 

0.2% 

0.3% 

Cspin(3:) 

= _o.066x) 

1-loop 7 "P“ 





Errors: 

The last four columns give the relative error for the parameterizations: 

• is the maximal deviation of the parameterization from the numerical data for the conversion functions 
at 3-loop precision, 

• is half the maximal difference between the numerical data at n and (n — l)-loop order, 

•^^6, MS = Apert, (m* = 4.25GeV), 

•^™C,MS = Apert, (m* = 1.2GeV). 


00 


Table 2.3: Parameterization of matching factors, x = 1/^(Mq/Aqcd) < 0.6. 
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the matching coefficients of the last section one obtains the following relations 
between meson decay constants in the full theory (l.h.s.) and in HQET (r.h.s.): 


IpS — Fps(mps)ymps 
Cps C'ps(Mq/Aqcd) 


tS'RGI + 0(1/Mg), 


Yy — Fv(mv)ymv 
Cv Cy{MQ/A qcd) 




(2.39) 

(2.40) 


In the static limit no interactions of the gluon held with the heavy quark’s spin 
survive and therefore 


/fvstat 

^PS, RGI 

/fvstat 

^V,RGI 


= 1 


(2.41) 


holds. which is the renormalization group invariant of the matrix el¬ 

ement dehned in table 2.2, has been computed non-perturbatively in [?] in the 
static approximation. 

For the ratios of decay constants one expects a behavior like 


R — Fpg(mpg) y'mps _i_ 

Cps/v Fv(mv) ymv C’pg/vCAfg/AQCD) 


1+0(1/Mg). 


( 2 . 42 ) 


Furthermore, HQET predicts the relation 

mx = mQ -I- A -I- ^ ^ Am^ -|- O^l/mq) for X = PS, V. (2.43) 
Q 


between the heavy-light meson mass mx and the heavy quark mass [?], where 


Am^ 


—Ai -|- 2 


J(J + 1 ) 



(2.44) 


J is the total spin of the meson, i.e. J = 0 for pseudo scalar mesons and J = 1 
for vector mesons. A is a parameter, that describes the properties of the light 
degrees of freedom in the background of the static color source provided by the 
heavy quark and and Ai oc (X| —Q{iD)'^Q\X) and A 2 oc The mass splitting 

gg:" c.,.7m;7aTo) = ( 2 . 45 ) 

where the lowest order contribution comes from A 2 , is therefore expected to van¬ 
ishes in the limit Mq —> cx). 

By producing data for the l.h.s. of these equations from relativistic lattice 
QCD for the range of heavy-light meson masses accessible to current lattice sim¬ 
ulations, it is the scope of this work to try to 

• estimate, down to which heavy quark mass in heavy-light meson systems 
observables scale with 1/Mq without sizeable contributions from the higher 
orders 
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• check the compatibility of the relativistic simulations with the static ap¬ 
proximation from an interpolation in 1/Mq 

• obtain a value for the decay constant and the mass splitting for a meson 
containing a 6-quark from an interpolation in the mass, i.e. including the 
prediction from HQET in the static limit. 

• estimate the order of the 1/Mq corrections to the static limit from a ht- 
ansatz a la 

Fps,v _ Qi 
Cps,v ° Mq 






Chapter 3 


Masses and meson decay 
constants on the lattice 


It has been shown in the previous chapter that the validity of HQET can be tested 
by exploring the mass dependence of meson observables in QCD. At their physical 
point such observables are also an important input for the phenomenology of the 
Standard Model. 

A particularly suitable framework in which a non-perturbative determination 
of mesonic observables is possible is the Euclidean QCD Schrddinger Functional 
on the lattice [?,?]• It has been demonstrated, that its Monte-Carlo simula¬ 
tion allows for the determination of mesonic observables with smaller statistical 
fluctuations than with conventional methods like lattice QCD on the torus. In 
addition, systematic errors introduced by excited states can be estimated more 
reliably [?]. 

First, the 0(a)-improved QCD Schrddinger Functional will be introduced in 
this chapter. Then, the meson mass and the decay constant of pseudo scalar 
mesons and vector mesons, and also the quark mass will be expressed in terms of 
renormalized and improved quark bilinear currents, whose expectation values can 
be evaluated in a Monte-Carlo simulation of the Schrddinger Functional. Finally, 
expressions of these currents in terms of quark propagators will be derived for 
the direct implementation in a computer program. 

3.1 The Schrodinger Functional - geometry and 
fields 

The QCD Schrddinger Functional on the lattice is the partition function 



( 3 , 1 ) 
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discretized on a hyper-cubic 4-dimensional Euclidean space-time cylinder 

Te = x/a e Z; 0 < Xo < T; 0 < Xfc < L; /c = 1, 2, 3} (3.2) 

with boundaries in the time direction. Here, is the QCD action dis¬ 

cretized on Te and will be specihed in the next section. The integration in (3.1) 
is the short form of 

f = [ \\d'ip{x)d'ip{x)dU{x). (3.3) 

J ^ 

The quark helds '4>a,A,a{x) are assignments of Grassman numbers to each lattice 
site x G Te and carry the flavor-, Dirac- and color-indices a, A, a respectively. 
The gauge helds U{x, p) G SU(3) are associated to the links between two adjacent 
sites {x,x + fi), fi being a unit vector in the /i-direction. 

In the spatial directions, the gauge helds U (x, /r) obey periodic boundary 
conditions 

U{x + Lk,k) = U{x,k), /c = l,2,3. (3.4) 

while the fermion helds 'ijj{x) are 6'-periodic 

^jJ{x + Lk) = e*®'='0(a^) and '0(x -|- Lk) = e“*®'='0(x), /c = 1, 2, 3. (3.5) 


The functional explicitly depends on Dirichlet conditions at the boundaries 
in the time direction. In particular, the gauge helds on the time slices Xq = 0 
and Xq = T are set to identity matrices^ 


G(x, fc)| 3 ,o=o = exp{Gfc} = U{x,k)i^^=T = expjC'fc} = Igxs iork = 1,2,3. (3.6) 

As the Dirac equation is a hrst order diherential equation, only two of the four 
components of the Dirac spinors 'ipix) on the boundaries have to be prescribed [?]. 
With = |(1 ± 7 o), the boundary conditions are 

P+ij{x)\^o=o = p{x), P_V'(x)po=T = p'{x), 


V'(x)P_Po =0 = p{x), '0(x)P+Po=T = p'(x). 


(3.7) 


3.2 Lattice action 


As in the continuum, a generic lattice action for QCD consists of a gauge action 
and a quark action S'f[P, ijj, ip], describing the dynamics of massive quarks 
coupled to the gauge held, 

_ S[U,ip,ip ] = Sg[U] + S^[U,iP,iP]. (3.8) 

^ Other choices of the boundary conditions have been used to define the running coupling 
constant in the Schrodinger Functional scheme, that scales with the size of the lattice volume 
like ^l=l/L [?,?]. 
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Wilson suggested the plaquette action, which in the version adapted for the 
Schrodinger Functional reads 

Sg[U] = - U{p)} . (3.9) 

9o „ 


The sum runs over all oriented plaquettes p. A plaquette is the smallest closed 
loop of link variables U{p) = U{x, p)U{x + fi,u)U~^{x + u, p)U~^{x,i'). The 
coefficient w{p) will be specihed in the next section. 

The quark action is 


{D + mo)'ll:(x). 

X 


Here, D is the Wilson Dirac operator 

0 = 5{v(V; + V,)-aV;V„)} 

with the covariant lattice forward and backward derivatives 

+ afi) - tplx)} , 


( 3 . 10 ) 


( 3 , 11 ) 


( 3 , 12 ) 


with Oo — 0 and —tt < d,, < tr (h — 1,2 ,3), is a phase factor that 
incorporates the spatial boundary conditions (3.5) for the fermion helds. Also 
the covariant lattice derivatives acting to the left side. 


= ^{Xl4i(x + a(i)U 

( 3 , 13 ) 

_ i - * _ _ 

- afi)U{x - afi,p)} , 

can be dehned. They will be useful at a later point. 

With the hopping parameter k = {8 + 2amo)~^, the fermion helds can be 
rescaled like ^l:{x) —> \/^ip(x) and 'ijj{x) —> \/^il){x). The fermionic part of the 
action can then be rewritten as 


Sf[U,%Ij,'iIj] = 'il!{x)M'il:{x), (3.14) 

x,y 


with 


M^{x) = ij^ix) - {Xf,U{x,p){l-'yf,)'ilj{x + afi) 

/ j .=0 

+ A*f/“^(x - a/i, ^)(l + 7^)V'(x - a/()} . 


( 3 . 16 ) 
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3.3 O(a)-improvement 

In Wilson’s original formnlation of lattice fermions [?], discretization effects dne 
to the finite lattice cut-off vr/a turn up at leading order in a. At the same time, in 
the Schrodinger Functional, also the Wilson gauge action is affected by 0(a) cut¬ 
off effects due to the presence of the boundary. Using an 0(a)-improved action 
and 0(a)-improved expressions for observables, lattice-discretization effects can 
be reduced (ideally removed) at leading order in a and the approach to the 
continuum limit is accelerated. 

The basic idea was first formulated by Symanzik [?,?], who showed for the 
^"‘-theory, that the corresponding lattice held theory can be described by a local 
effective continuum held theory with the action 



Co is the generic continuum Lagrangian. The terms Ck {k = 1,2,3,...) are 
combinations of local operators of dimension 4 + k, all of them respecting the 
exact discrete symmetries of the lattice theory. In the ehective continuum theory, 
the lattice helds are represented by the renormalized ehective helds 


(l)es{x) = 0o(x) -F a0i(x) a^(t) 2 {x) + .... (3.17) 

The helds 0fc(x) (fc = 0,1, 2, 3,...) also must have the appropriate dimension and 
respect the same lattice symmetries. 

The 0(a)-improved theory can now be obtained by adding suitable counter 
terms with properly tuned improvement coefficients to the lattice action and to 
the lattice helds, such that the 0(a)-terms in the ehective continuum action and 
the ehective continuum helds are canceled. In the cases of interest for this work, 
only observables on the mass shell will be considered. Thus, on-shell improve¬ 
ment can be applied, where the number of counter-terms, that are necessary 
to improve the action and the helds at 0(a) can be reduced by exploiting the 
classical equations of motion. 

0 (a)-improvement has shown to be an efficient tool for lattice-held theory [?]. 
The improvement coefficients have been determined non-perturbatively [?,?] in 
most of the cases, or otherwise are available from perturbation theory [?,?,?]. 

One thing to mention here is, that Symanzik’s proof, that 0(a)-improvement 
works, relies on perturbation theory. On the lattice one therefore always assumes, 
that 0(a)-improvement is applicable beyond perturbation theory. Experience 
from lattice simulations of QCD however support, that this assumption is justihed 
[?,?,?,?,?]. 
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3.3.1 Improved action 

To cancel the 0(a)-effects in the Wilson SU(3) gauge action in the Schrodinger 
Functional, one has to choose the weight w{p) in (3.9) as [?,?] 

{ |cs(fi'o) if P is a spatial plaquette at Xq = 0 or xq = T, 

Cf((7o) if P is a time-like plaquette attached to a boundary plane, 

1 elsewhere. 

(3.18) 

c<j(po) can be dropped for simulations with vanishing background held, as it will be 
the case throughout this thesis, while Ct{go) has been determined in perturbation 
theory at 2-loop order [?] and is given in appendix C. 

For the 0(a)-improvement of the fermionic part of the QCD Schrodinger 
Functional, volume and boundary counter terms 6Sy[U,'ijj,'ip] and 6Sb[U,'ijj,'ip] 
have to be added to the action which then reads 

V', tP] = Sf[U, 'll;, ip] + 6Sb[U, ^P] + 6Sy[U, Pj]. (3.19) 


The superscript I from now on indicates 0(a)-improvement. The boundary im¬ 
provement term is 


5Sb[U,'ip,'ip] = +0',{x) 

(q(po) -i)\dt{x)- d((x) 


(3.20) 


where 


Os{x) 

= P(^)^7fc(Vfc +Vfc)p(x), 

(3.21) 


= p'(^)^7fc(Vfc +Vfc)p'(f), 

(3.22) 

dt{x) 

= lpj{x)P+'\/lp>{x) + ip{x) Wo P-P’ix)^ and 

L J \xo=a 

(3.23) 

o[{x) 

= \'ip{x)P^Wo'ip{x)+'ip{x) Wo P+'ipix)} 

L J \xo=T-a 

(3.24) 



(3.25) 


The volume- or Sheikholeslami-Wohlert-term [?], has the form 

T-a 

5Sy[U,^P,^P] = csw(fi'o)^cr^i/ ■ F^y{x)ip{x). (3.26) 

XQ=a X 


= f b/.!) Ii] and is the symmetric lattice held strength tensor 


(3.27) 
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with 


Qfiuix) = U{x, ^)U{x + afi,i')U{x + aO, 

+U{x, i')U{x — aft + ai>, fi)~^U{x — afi, h')~^U{x — afi, /i) 

+U{x — afi, ^)~^U{x — afi — aP, h')~^U{x — afi — aO,fi)U{x — aP, v) 
+U{x — aP, h')~^U{x — aP,fi)U{x + afi — aP, h')U{x, 

(3.28) 

Qfj,u{x) is the sum over four plaquettes in a hyper plane, all having the site x in 
common and therefore resembling a clover leaf. 

The improved fermion action can equivalently be expressed in terms of an 
action with the improved Dirac operator 

= D + dD^ + 6Db, (3.29) 

where 

6D^'if{x) = acsw{9o)^(^fiu ■ F^,,{x)'ijj{x) (3.30) 

and 


SObfiix) 


(ci(^o) - 1) 


1 

a 



[lf{x) -U{X- 


aO) ^P+ifi^x — aO)] 


(3.31) 


+ ^xo=T-a - U{x, 0)P_'lf{x + aO)] 

As for the unimproved Dirac operator, a parameterization in terms of the hopping 
parameter k is possible, namely 


+ mo) ^|J{x) 


1 

2k 


{M + SM)if{x) 


^M^if{x) 

2k 


(3.32) 


with 6M = 2k {SDy + 6Db). 

The improvement constant Ct{go) has been determined in perturbation theory 
at 1-loop order [?] and csw(5'o) has been determined non-perturbatively [?]. Both 
constants are given in appendix C. 


3.3.2 Improved fields 

As it will be explained in section 3.6, all phenomenological observables that have 
been determined in this work, can be expressed in terms of quark bilinear currents. 
Of particular interest are the iso-vector axial current Af{x) = '^(a;) 7 ^ 75 ^- 0(^)5 
the iso-vector vector current Vff{x) = ip{x)'y^^if{x) and the iso-vector pseudo 



3.4. FERMIONIC OBSERVABLES IN THE SCHRODINGER 
FUNCTIONAL 
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scalar density P“(a;) = Except for the pseudo scalar density, these 

currents get contributions from dimension hve operators under improvement [?]. 
The improved axial current is 


= A“(a:) + CA{go)adf,P'^{x), 

(3.33) 

where 



(3.34) 

is the symmetric lattice derivative. For the vector current V))*(x) 
one adds the symmetric derivative of the tensor current 

= y{x)'yi,^ip{x), 


(3.35) 

{yyYix) = v;^{x) + ^cy(c/o)a4T;^ 

(3.36) 


ca(5'o) and cv{go) have been determined non-perturbatively in [?] and [?,?], and 
their parameterizations in terms of the bare coupling^ are given in appendix C. 


3.4 Fermionic observables in the Schrodinger 
Functional 

Summarizing the last section, the 0(a)-improved Euclidean QCD Schrodinger 
Functional with the lattice action 

S^[U,^P,^P] = S^^[U] + Si[U,^l^,^P]. (3.37) 

and vanishing background field is 

Z[p,p,p',p']= [_ (3.38) 

Sq and Sp are the improved Wilson gauge and fermion action introduced before. 
A source term 

= Y{vix)Hx)+'ilj{x)vix)} (3.39) 

0<3;o<T X 

will now be added to the action. One can then identify 


^{x) 

s 

= 

S 

5fj{x) ’ 

Sr/(x) ’ 

C{x) 

5 

C(^) = 

s 

5p{x )’ 

Sp{x )’ 

C(x) 

S 

C'(^) = 

<5 

5p'{x) ’ 

5p' (x) 


^In the case of cvigo)-, the parameterization derived in [?] has been used here 
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where C(^)) C(^)) ('(^) ('(x) can be interpreted as boundary fields. Ex¬ 

pectation values of operators that are polynomials O in these fields are defined 
as 


(O) 


{14 


4,1/ 


(9 Q-'S^lU,i/j,tp]-Sslip,ip,V,v] 


lp=p=p'=p'=r/=i)=0 


( 3 . 41 ) 


With O being the polynomial O, where all fields replaced by the associated 
functional derivative, the expectation value can be obtained as 


(O) = ld\YiZ[p,p,p\p\ri,r]\\ . (3.42) 

l. J \p=...='q=0 

For the following, only the fermionic contribution to {O) on a given gauge 
background, 

[C)]p[f/] = ^ (3.43) 

with 

Zf[D,p,p,p',p',4,4] = / ^-sUuZM-Ss[M,v,v] (3.44) 

will be considered. The relation to the full expectation value is given by 


{O) - {[(!1]f)g> 


( 3 . 46 ) 


where the subscript G indicates, that the expectation value has to be evaluated 
with respect to the path integral of the quenched theory. 


3.5 Renormalization 


Usually, a mass independent renormalization scheme is employed [?] in which the 
renormalization constants are independent of the quark mass and therefore the 
renormalization group equations, which describe the scale dependence of renor¬ 
malized quantities, take a particularly simple form. In this scheme, the theory is 
parameterized around a critical line Kcritido) for which the subtracted bare quark 
mass 


arrig 


^crit (40 


(3.46) 


and also the renormalized quark mass is zero. 

All the bare parameters and fields receive a multiplicative renormalization 
factor. However, when taking the continuum limit, uncanceled 0(amq)-cutoff 
effects arise [?]. In the improved mass independent renormalization scheme, one 
subtracts these effects by a multiplicative counter term of the form 1 -|- b{gQ)amq 
and defines [?] the improved renormalized subtracted quark mass (with flavor 
index i = s,c,b) 


= Zm{go)mq^i = Zm{go)mq^i{l + bm{gl)amq^i). 


(3.47) 
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The renormalization of the heavy-light quark currents can be discussed con¬ 
veniently after introducing the off-diagonal iso-matrices 

r^ = r^±ir^. (3.48) 

Then, the improved renormalized axial current 

(^J)R(a^) = ^a(^o)(1 + bA{go)^{amg^h + amg^i)){A^y {x), (3.49) 

the vector current 

= Zv{go)y + bv{go)^{amg^h + amgy){Vyy {x), (3.50) 

and the renormalized pseudo scalar density 

= Zp{go,fi){l + bp{go)^{amg^h + amgy)P^{x) (3.51) 

can be dehned, where h and I indicate a heavy and light flavor respectively. The 
mass renormalization constant Zm{go) and the scale dependent renormalization 
constant Zp{go,ix) will only be needed in terms of the composite renormaliza¬ 
tion factor Z^go) = Zm{go)ZA{go)Zp^{go) (cf. section 3.6.3). Z{gQ) has been 
determined non-perturbatively in [?] and ZA{go) and Zyi^go) in [?] and [?]. Their 
parameterizations in terms of the bare coupling constant go have been summa¬ 
rized in table C.l in appendix C. The improvement constants bA^go), by {go), 
bp{go) and bm{go) have to be tuned in order to compensate for the 0{amg)- 
terms. 5/i(5'o) and by {go) have been determined non-perturbatively in [?] and 
bm{go) ia [?]• bp{go) will only be needed in the difference bA{go) — bp{go), which 
also has been determined in [?]. All the 5-factors are given in appendix C in table 
C.2. 

3.6 Meson masses, decay constants and quark 
masses in the Schrodinger Functional 

In this section, expressions for the pseudo scalar (X=PS) and the vector me¬ 
son (X=V) decay constants Ex and the corresponding meson masses mx will be 
derived. In the the Schrodinger Functional, this can be done in terms of corre¬ 
lation functions of quark-bilinear currents O carrying the appropriate quantum 
numbers, 

X = PS : O = Aq{x) = y{x)'^ol5^'4^{x) and 


X = V : O 


vy{x) = %i){x)-^i'^iij{x). 


(3.52) 
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Furthermore, expressions for renormalized quark masses will be given in this 
framework. 

The derivation is based on the transfer matrix formalism which has been 
formulated for the QCD Schrodinger Functional in [?]. Although the derivation 
does not carry over rigorously to the improved theory, universality implies that 
the construction is valid there as well and that one can replace the currents (3.52) 
by the improved currents (3.49) and (3.50) [?]. 

3.6.1 Correlation functions in the Schrodinger Functional 

The lattice Schrodinger Functional for QCD (3.1) [?] can be represented as the 
QCD transition matrix element^ 

Z[C,C\p,prp\p'] = (*'| (e-®)^PK) (3.53) 

between an initial state |i) and a final state {i'\ (at times Xq = 0 and Xq = 
T respectively). The exponential with the QCD-Hamiltonian H, can be 

interpreted as the transfer matrix connecting two adjacent time slices. P is a 
projector constraining the dynamics to the gauge invariant subspace of the whole 
Hilbert space. 

Relevant for this work are the eigenstates |n, q) of the Hamiltonian which are, 
next to the energy quantum number n, specified by the set of quantum numbers 
q = {J,C,P,mh,mi) (total angular momentum, parity, mass of the heavy and 
light quark respectively). 

The pseudo scalar Dg-meson for example can be specihed by the quantum 
numbers q = (0, ±, —,mc,ms), the one of the corresponding vector meson D* by 

q = (1, ±, 

The \n, q) build a complete set of energy eigenstates of the mesonic sector of 
the Hamiltonian IH with 


\n,q), n = 0,1,2,..., (3.54) 

M\n, q) = E^\n, q) (3.55) 

and the normalization 

(n', q'\n, q) = 5n,n'Sq,q'. (3.56) 

In the case of a Schrodinger Functional with vanishing fermion and gluon 
boundary fields, 

K') = K) = Ko) (3.57) 

holds, with |io) carrying the quantum numbers of the vacuum. 

^The boundary fields C and C are mentioned for completeness here and will again be set 
to 0 in the following. 

'’’As electromagnetic effects only play a minor role here (cf. [?]), the associated charge will 
be neglected. 
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C' r'o C' 


Xo = T 



Figure 3.1: The Schrodinger Functional correlation functions fo and Jq. 


In the following, the focus will be on two types of matrix elements. The hrst 
matrix element is 

foixo) = -Z-'i(*o|e-(^-"°)®PO(T)e-"°®P|zx). (3.58) 

It measures the correlation between the operator 0(T) at time Xq with the initial 
and hnal states |ix) and (^ol- The physical picture is that of a meson with 
quantum numbers gx being created from fermionic boundary sources at time 
Xo = 0 and being annihilated at some later time Xq by the operator 0(x) (cf. 
hgure 3.1). The second matrix element is^ 

/o = (3.59) 

It measures the correlation between the initial state |ix) and the hnal state 
after time evolution over the time extent T and thereby represents a meson state, 
traveling through the space-time from the boundary at xq = 0 to the boundary 
at xo = T. In practice, the states |zx) and (zxl are created by dimensionless 
mesonic boundary operators in the zero-momentum projection. 


O 


a _ a® 

O — 


Ecwnyc®i«.o. 

^,y 


o 


a! _ a® 

O — 




xo=T- 


(3.60) 


and the operator 0(T) typically represents a bilinear quark current of the form 

0“(a;) = '0(a;)roy^/^(x). (3.61) 

The Dirac matrices Tq are collected in table 3.1 for the axial current (A), 
the vector current (U), the pseudo scalar density (P) and the tensor current (T). 
They determine the desired quantum numbers of angular momentum and parity. 

^In the literature on the Schrodinger Functional, these matrix elements are also often refered 
to as fi = fj, and fci = /^. 
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0 

A 

V 

P 

T 

To 

7o75 

li 

75 

o'ij 

V' 

0 

75 

-\li 

75 

-bi 


Table 3.1: Dirac structure for the currents 0{x) = f[ 7 i, 7 i])- 


The matrix T^ = 7075 for example corresponds to {J,P) = (0,—), while the 
matrix Ty = :^ 7 ^ corresponds to {J,P) = (1, —). The SU( 2 ) flavor structure of 
the current is given by the matrices (a = 1, 2, 3, cf. appendix A). 

The matrix elements fo{xo) and /q can then be identified with the expecta¬ 
tion values 

foixo) = -y (0+(xo)(Po) for 0 = A,P,V,T (3.62) 

and 

fT = foTX = A,P (3.63) 

where again the flavor off-diagonal fields (cf. (3.48)) have been used. In order 
to arrive at expressions for the meson mass mx and the meson decay constant 
Ex in terms of these correlation functions one first inserts twice the identity 
1 = Z] fofo (3.58). This yields 


foixo) =—Z (*o|e g)(n, g|0|n', g')(n', g'|e ^°’®P|ix) 

q, q' 

n, n' 

~ _^-iY|y|g-(r-xo)ep(|o^O)(0,0| + |0,gx)(0,gx| 

+ |1,0)(1,0| + |l,gx)(l,gx|) O (|0,0)(0,0| + |0,gx)(0,gx| 

+ |l,0)(l,0| + |l,gx)(l,gx|) e-"oEp|*x)}. 

(3.64) 

In the second step, all but the ground state, the first excited state = 0,1 
and q = 0 , gx have been neglected, since their contribution will be suppressed 
exponentially. With 


Z^(*o|P|0,0)(0,0|P|*o)e-''™'^ 


(3.65) 
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one arrives at 
fo{xo) ^ 


X 


1 pEf^T 

2 (io|P|0,0)(0,0|P|io) 

I (*o|P| 0 , 0 )( 0 , 0 | 0 | 0 , gx)( 0 , 


+ (*o|P|0, 0)(0, 0|O|l, gx) (1, gx|PKx)e-''™U-o)e-®?^-o 

+ (*o|P|l,0)(l,0|0|0,gx)(0,gx|PKx)e-<’(^-"°)e-®o^"° 


+ (zo|P|l,0)(l,0|0|l,gx)(l,gx|PKx,)e-'''°'(^-"°)e-®i^"°|. 


(3.66) 


In order to write this in a more compact form, the following abbreviations can be 
introduced. First, rug^ = — Eq and rriQ = E^ — Eq are the mass of the ground 

state meson and the glueball respectively, A = Ef" — Eq^ is the gap energy in 
the meson sector. Furthermore one dehnes the matrix elements 


p = 

(0, gx P 7x) 

(0,0|P|7o) ’ 

(3.67) 

Px = 

(0,0|O|l,gx)(l,gx|PKx) , 

(0,0|0|0,gx)(0,gx|P|7x) '''' 

(3.68) 

pI = 

(7o|P|l,0)(l,0|0|0,gx) 

( 7 o|P| 0 , 0 )( 0 , 0 | 0 | 0 ,gx)' 

(3.69) 



(3.70) 


fo then takes the simple form [?] 

foi^o) ^ (0, 0|0|0, gx) , 

(3.71) 

which decays with the mass of the meson. For Jq one derives 

rT 1 l(^x|P|0,gx)P T{El^-E°g) _ 1 2 (o 72 ') 

2 |(7o|P|0,0)|2 2^ ^ ’ 


The correlation function fo{xo) contains contributions from higher excited states 
for small Xq, which decay exponentially with the gap energy A. For large Xq, 
bound states of gluons, the so called glueballs contribute exponentially enhanced 
with their mass rriQ. 

For the numerical implementation it is very convenient to also dehne the 
backward correlation functions 

go{T-xo) = f(0+'0-(xo)) forO = A, U, 


go{T -xq) 


f(0+'0-(xo)) for 0 = P,T. 


(3.73) 
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(next to the forward correlation functions fo)- They are constructed in the same 
way as fo{xo) but with the meson sources at the time slice at Xq = T and serve 
the same purpose, provided that the background held is zero. In a Monte-Carlo 
simulation, both the forward correlation functions and the backward correlation 
functions will be evaluated on the same sample of gauge conhgurations. Although 
the resulting data will be correlated, the statistics can be increased in this way. 

The improvement of the axial current (3.49) and the vector current (3.50) is 
taken over from [?] to dehne the improved correlation functions 

fA{xo) fA{xo) = fA{xo) + acAigo)dofp{xo) (3.74) 

and 

fv{xo) fvi^o) = fvixo) + ^acv{go)BfT{xo) (3.75) 

and analogously for the correlation functions go(T — xq). 


3.6.2 The meson mass and the meson decay constant 

The meson mass mx can be extracted as the effective mass. 


meff(Xo f) 


iln 


\fo ho+a.) ) 


2sinh(aA/2) qx -xpA _ 2 sinh(amG/2) Q (T-xn)mr, \ 

— TJy^e J. 


The connection to the meson decay constant, which is dehned as [?] 

Zo(0, 0|O|0, gx) = Fxmx(2mxT^)“^'^^ 


(3.76) 

(3.77) 


is then given by 


X 


-2Zo(mxT3)-V2 g(xo-T/2)mx/^ 

V fo 


(3.78) 


Here, Zq is the renormalization constant of the current O and (2mxT^)“^/^ in 
(3.77) is the conventional normalization of one-particle states. 

In terms of the backward correlation functions, the effective mass and the 
decay constant can be expressed by replacing /^(xo) with gpiT — xq) in (3.76) 
and (3.78). 


3.6.3 Quark masses 

On the lattice, there exist various ways to dehne quark masses, which differ at 
hnite lattice spacing but converge in the continuum limit. One dehnition has 
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already been introduced in (3.47), where the renormalized valence quark mass is 
given in terms of the hopping parameter. Another way to obtain the quark mass 
is guided by the PCAC relation in the continuum 


= 2mP“, (3.79) 

where m is the quark mass. The average bare current quark mass in the Schrodinger 
Functional can be written as [?] 


mhi = 


-((9o + do)fA{xo) + CAadQdofp{xo) 


/fp{xo). 


(3.80) 


The sum of the renormalized valence quark masses can then be obtained by 
multiplicative renormalization of the axial current (3.49) and the pseudo scalar 
density (3.51) [?], 


rriR^h + mYi,i 


^ ZAil+bAi9o)^iamq^h+arnq^l)} , 

Zp(l+bp (go) I (amg^h+arrig^i )) 1 


= 2ZaZp\1 + (bA(go) - bp(go))^(amg^h + amq^i))mhi + 0{a?). 

(3.81) 


3.6.4 Correlation functions, propagators and symmetries 

For their implementation in a Monte-Carlo simulation program, the correlation 
functions fo{xo), goiT — Xq) and /q must be expressed in terms of contractions 
of quark propagators which will be derived in the following. 

Writing down the flavor indices explicitly fo{xo) reads 


fo{xo) = 




E {'ilJh{x)ToMx) Ciiy)'^'oCh{^)- 


y,z 


(3.82) 


for O = A, P, V, T. Therefore one can write fo as 

fo{xo) = ^E(Tr{ [Ch(^)^/i(a;)]FTo [V^z(x)CKy)]Fro}>G’ (3.83) 


y,z 


where the trace is over Dirac and color indices. Analogously one obtains 

goiT-xo) =±^^^{TT{[ijh{x)Ch{y)]p'^o[Cii^Mx)]F'^o})Q, ( 3 . 34 ) 


y,^ 


where the sign has to be chosen according to (3.73), and 


fo = 




(Tr{[a(^ia(^)]Fro [c;(^)CK^)]Fro}>G- 


v,w^y,z 


2 


(3.85) 
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Applying the functional derivatives (3.40) to the fermionic generating functional 
In Zf[D, p, p, p', p', f), r;], the propagators [ ■ jp in (3.83), (3.84) and (3.85) can now 
be written down explicitly: 




(3.40) 


CtP-U{z,Q)S{z,x)\^^=a, 


|p=..-r,=0 


(3.86) 


[V'(a:)C(p)]p 

= CtS{x,y)U-^{y,0)P+\y^=a, 

[V'(a:)C'(p)]p 

= CtS{x,y)U{y,0)P-iy,=T-a, 

[C'(^^(x)]p 

= CtP+U-^{z,0)S{z,x)l^a=T-a, 

[C(50C'(n))]p 

= c‘fP-U{z,0)S{z,v)U{v,0)P_\^Q=a,vo=T-a, 

[C'(hi)C(y1)]p 

= cfP+U-^iw, 0)S{w, y)U-^{y, 0)P+\^^=T-a,yo 


(3.87) 


Inserting these expressions into /o(a^o) one gets, denoting the quark flavor by the 
subscripts h and I respectively. 


foixo) = {P-U{z,0)S,ix,z)roSi{x,y)W{y,0)P+To}) 




G| 




E (Tr I P+U{z, Q)Si{x, z)-i,T'oSi{x, y)U\y, 0)P+ro75} 

y ,2 ' 


G| 


!/0=^0 = 


= n (^sl(xy„r'oS,{x)ro'r,))^, . 

(3.88) 

Here, the summation over the boundary fields has been included into the propa¬ 
gator 

^(x) = Cia^^^(a:,p)P"^(p,0)P+|^^^^, (3.89) 

y 

which is the propagator of a zero-momentum quark state at po = 0 to a space- 
time point X in the interior of the Schrodinger Functional. In the same way, 
go{T - xo) is 


go{T-xo) = (Tt (^Rl{x)j5T'oRi{x)To^5^^^. (3.90) 

with 

R{x) = Cta^^S'(x,p)P(p,0)P_|yo=r_a, (3.91) 
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being the propagator of a zero-momentum quark state aXy^ = T to the space-time 
point X. The boundary-to-boundary correlation function Jq is 


s2 6 7-3 , , ^ >1 \ 

= 5](iv{st(„)[/(i7,o)p+75r'oP+c/-‘(«>,o)s,(«>)ro})^, 

Z ^ \ ^ ) ! G|vo=T—a, wo=a 

v,w 

(3.92) 

Introducing the boundary-to-boundary quark propagator 

St = CtY,P+U-\x,Q)S{x)\,,=T-a (3.93) 


it takes the form 


Observing that 




S{x)P. = R{x)P+ = 0 

in the chiral basis for the gamma matrices, the relations 


(3.94) 

(3.95) 


S'ai(x) S_az{.x) = 0, Sa2{,x) + S_A4{x) = 0, 

Rai{x) - Ra3{x) = 0, Ra2{,x) - Raa{x) = 0, 


hold, where A is a Dirac index. Thus, all correlation functions introduced so far 
are completely determined if only half of the Dirac-components of the forward- 
and backward propagators S{x) and R{x) are known. 

In the numerical simulation, these components of S{x) and R{x) have to be 
determined on each gauge background of the Monte-Carlo history. The equations 


[D^ + m)S{x) = + ^) S{x,y)U ^{y,Q)P+\y^=a 

y 

= cta~^5^Q^aU~^{x,R)P+ 


(3.97) 


and 


+ m) R{x) = Cttt ^6xQ,aU{x,0)P^, (3.98) 

or equivalently 

MG(x) = 2KCta-^6,,,aU-^ix,0)P+, 


(3.99) 


JVRR{x) = 2KCta ^6xo,aU{x,0)P- 

therefore have to be solved for S{x) and R{x), using a matrix inversion algorithm, 
the stabilized biconjugate gradient [?] for example. 




Chapter 4 
The PC-Code 


For a major part of the simulation parameters, the production runs could be ac¬ 
complished on the APEMille supercomputers at the John von Neumann Institute 
for Computing [?]. In particular, the simulations for the scaling study with the 
bare couplings (5 = Q/Qq = 6.0, 6.1, 6.2, 6.45 ^ were all done within this framework 
using the ALPHA-collaboration’s approved program code for the simulation of 
quenched QCD in the Schrodinger Functional. In simulations for the lattice with 
the hnest resolution, i.e. at /3 = 6.7859, the number of lattice points L/a and 
therefore the amount of memory space that needs to be allocated to store the 
associated helds gets very large and the available memory on the APE-computer 
is not sufficient any more. 

An alternative supercomputer with the appropriate specihcations was the 
newly installed system of 26 IBM pSeries 690 eservers at the Norddeutscher 
Verbund fur Hoch- und Hdchstleistungsrechnen (HLRN). Each of these servers 
has 32 IBM-Power4 processors (with a peak performance of 5.2 GFlop/s each) 
that share between 64 to 256 GByte of memory. 

As the lattice-gauge code on the APE-computer is written in a proprietary 
language (TAG), a new simulation code had to be obtained. To this end, the 
lattice code by the MIMD^ Lattice Computation (MILC) Collaboration [?], pub¬ 
lished under the GNU General Public License [?], seemed to be an appropriate 
choice. It consists of a set of routines written in C for doing simulations of four di¬ 
mensional SU(3) lattice gauge theory. The version which the collaboration shares 
on the Web, already comprises many of the features that are necessary for the 
realization of this project. Among those are 

• Schrodinger Functional boundary conditions, 

• Wilson pure gauge code, 

^Here, denotes the bare gauge coupling. 

^The aforementioned IBM-computer has a MIMD architecture {Multiple Instructions, Mul¬ 
tiple Data), as opposed to a computer architecture like the one of the APEMille which is a 
SIMD {Single Instruction, Multiple Data) computer. 
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• Wilson Dirac operator, 

• 0 (a)-improvement, 

• stabilized biconjugate gradient inverter, 

• platform-independence, 

• MPI-based parallelism (MIMD). 

This chapter discusses the MILC Code, all changes and improvements added 
to it and the extensive testing previous to the data production. Other topics are 
the resource allocation for the production runs and the corresponding applications 
for CPU-time at the HLRN. 


4.1 The MILC Code 

The two building blocks in the lattice simulation of gauge theories are the gen¬ 
eration of a Monte-Carlo series of gauge-configurations and the subsequent eval¬ 
uation of observables on the given gauge background. 

The gauge updates in the MILC Code are done with a quasi-heatbath gauge 
update by three SU(2) subgroups a la Kennedy-Pendleton [?] and Cabibbo- 
Marinari [?]. Also micro canonical over relaxation steps by doing SU(2) gauge 
hits have been implemented. The pseudo random number generator is the ex- 
clusive-OR of a 127 bit feedback shift register and a 32 bit integer congruence 
generator. It runs with a different seed on each of the parallel processors to avoid 
correlations. 

Based on the even-odd preconditioned Dirac operator [?], the quark propaga¬ 
tors on a given gauge background are evaluated using the stabilized biconjugate 
gradient algorithm [?]. 

Up to the random number generator, the implementation of the algorithms 
resemble very much the ones in the ALPHA Collaboration’s lattice code. How¬ 
ever, some peculiarities should be mentioned here. The MILC code uses the 
Weyl-basis (75 diagonal) Dirac matrices. In contrast to the standard convention, 
the projectors P± = |(1 ± 70 ) are implemented with opposite sign, that means, 
the projectors in the MILC code 

pMILC ^ 

The MILC code uses MPI-based (MPI stands for Message Passing Interface [?] 
) parallelization. The code has been designed in such a way, that the user in most 
cases does not have to bother about the parallelization, e.g. when implement¬ 
ing a set of correlation functions. Only at a deeper level, for example for the 
implementation of a new action, these issues have to be considered. Assuming, 
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that the number of lattice sites is divisible by the number of nodes, which is a 
power of two, the lattice volume is divided by factors of two in any of the four 
directions. Dividing the direction of the largest lattice extent is favored in order 
to keep the area of the surface minimal. Similarly, dividing directions, which 
have already been divided is preferred, thereby keeping the number of off-node 
directions minimal. 

The root project directory has the following listing: 

generic/ 

generic.clover/ 

generic.pg/ 

generic.schroed/ 

generic.wilson/ 

include/ 

libraries/ 

schroed.pg/ 

f_A/ 

The two main directories are schroed.pg/ and f_A/. They contain the main 
program hie for the gauge-update and the code for the evaluation of observables 
on a given gauge background. 

The generic/ directory contains high level routines that are more or less 
independent of the physics. Examples of generic code are the communication 
routines (MPI), random number routines, routines to evaluate the plaquette, etc. 
A set of slightly more specihc directories is generic.clover/, generic.wilson/, 
and generic.schroed/. They contain the implementation of the Dirac operator, 
the clover term, the matrix inversion routines or routines for setting up the lattice 
topology. All applications share the include/ directory, containing most of the 
header hies and the libraries/ directory containing low-level routines, mostly 
linear algebra routines. 

Detailed instructions on how to use the code can be found in appendix E. 
Further information about the MILC code in general can be found on the project 
web site [?] 

4.2 Changes 

The MILC code as it is available on the web had to be customized in a couple of 
places in order to conform to the requirements demanded by the project. 


4.2.1 The main program 

Next to setting up the lattice including the fermionic sources, initializing MPI- 
based inter-node communication and loading the gauge conhguration, the main 
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program control_cl. c contains a nested loop in which the inverter is called for 
a given set of indices like 

forCcolor = 0; color < 3; color++){ 
forCspin = 0; spin < 2; spin++)-[ 

for (kappa = 0; kappa < nuni_kap; k++)-[ 

*** CALL BiCG (color, spin, kappa)*** 

} /* kkappa */ 

for(kkappa = 0; kkappa < num_kap; k++){ 
fordkappa = k; Ikappa < num_kap; !++){ 

*** EVAL CORRELATION FUNCTIONS (color, spin, kappa) *** 

}- /* Ikappa */ 

} /* kkappa */ 

} /* spin */ 

} /* color */ 

In this loop-structure, the symmetries of the propagators pointed out at the end 
of section 3.6.4 have already been taken into account. The loop over the spin 
components only runs from 1 to 2. Only one color component of the propagator 
has to be stored in the memory at a time. Whether the code evaluates the 
forward- or backward-propagators has to be decided at compilation time and 
just affects, the way the fermion sources are set up at the boundaries. 

The contraction of the propagators that give the correlation functions given in 
3.6.4 have been implemented in the routines f_A. c, f_P. c, f_1. c, f_V. c, k_T. c, 
k_l. c. 

4.2.2 Double precision arithmetics 

The official MILC code with Schrodinger Functional boundary conditions is based 
on single precision arithmetics^. For simulations of large lattices, this may become 
a source of concern, especially when evaluating sums over fluctuating numbers 
over the whole lattice. To be on the safe side, double precision arithmetics have 
been implemented throughout the code. 

A compiler flag now allows for changing back to single precision arithmetics at 
compilation time (dehned in include/config.h). Thus, allowing at the price of 
reduced precision, to allocate only roughly half the memory necessary for double 
precision arithmetics. 

A better number precision results in a better precision of the inversion al¬ 
gorithm. This has been investigated with the following test. First, the Dirac 
operator applied to a quark source 0in was inverted on a gauge background 

^According to the IEEE-standard, a single precision number consists of a sign bit, eight 
exponent bits and 23 bits for the mantissa, while a double precision number consists of one sign 
bit, 11 exponent bits and 52 bits for the mantissa. 
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Figure 4.1: Exact residuum ||0in — 0out||/||0in|| plotted against the iterated 
residuum using single precision arithmetics (crosses) and the newly implemented 
double precision arithmetics (circles). The dashed lines indicate a residuum of 
10“® and IQ-i®. 


(L/a = 48® X 96, /3 = 6.7859, k = 0.117625). Then, the Dirac operator was 
again applied to the result, i.e., 

0out = dT(M“^0in)BiCGstab- (4.2) 


The inversion algorithm monitors the convergence with an iterated solver 
residual e [?]. After every iteration, it can be constructed from the change of the 
solution vector with respect to the previous iteration. Compared to the exact 
residuum 


Il0in 


0out 


Il0in 


( 03:.in */*3;,out) 

X 




^x,in 

X 


(4.3) 


one thereby saves one time consuming application of the Dirac operator in each 
iteration step. 

The solution 0out is accepted as soon as the iterated residual is smaller than 
the stopping criterion e. This computation has been repeated for a number of 
residuals e = 10“^,10“®,... ,10“^®. The expectation is, that for single precision 
arithmetics, a residual of ~ 10“® can be reached and ~ 10“®® in the case of double 
precision. 

The results shown in hgure 4.1 conhrm, that the implementation works in the 
desired way: The exact residual decreases with e and then saturates at about the 
corresponding arithmetic precision. 
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4.2.3 Performance 

A well-known disadvantage of the MILC code is the way the allocation of memory 
is organized [?]. All variables associated to a lattice site are stored in a local 
structure called site - metaphorically speaking this is a container for indices, 
gauge helds, auxiliary helds etc., dehned at each lattice site. It has the form 

typedef struct { 

/* coordinates of this site */ 
short x,y,z,t; 

/* is it even or odd? */ 
char parity; 

/* index of the site in the lattice array */ 
int index; 

/* Physical fields, application dependent, 
add or delete whatever is needed.*/ 

/* gauge field */ 
su3_matrix link[4]; 

/* spatial boundary links */ 
su3_matrix boundary[3]; 

/* temporary link variable for Field Major /* 
su3_matrix link_tmp[4]; 

/* Wilson complex vectors */ 

wilson_vector psi; /* solution vector */ 

wilson_vector chi; /* source vector */ 


} site; 

(c.f. lattice.h in the project directory). 

The sites on each node of the parallel machine are the elements of a local 
array which occupies a linear space in the node’s memory. One often refers to 
this type of memory arrangement as site major because one hrst has to point to 
a particular site before accessing the appendant variables. 

On the one hand site major is concise and user-friendly as it allows to easily 
add new variables, auxiliary helds e.g., that may become necessary when creat¬ 
ing new projects. On the other hand, it is very ineffective with respect to the 
communication between the main memory and the computer’s cache. 

One reason for this is, that during a lexicographic sweep through the lattice, 
one often asks for the value of a certain variable at the current site. As detailed 
before, not each variable on its own, but the array of sites occupies a connected 
address space in the memory. In each single step of the sweep, the address pointer 
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Figure 4.2: Sketch of site major and field major cache processing 


therefore has to be translated by the length of one site in address space, causing 
sizeable latency. 

Another reason is intimately connected with hardware prefetching, a common 
feature of modern processor architectures (e.g. IBM Power4), where the processor 
guesses which address space has to be prefetched into the cache for processing 
during the next clock steps. This feature only works properly if the variables to 
be processed in a loop are allocated in a reasonably homogeneous way. 

A third point is, that some processors allow for multiple data streams (e.g. 
IBM p690 up to 8 streams). That means that the processor can handle a number 
of simultaneous data interchanges - streams - between the main memory and 
the cache. The processor can identify streams only, if the data is contained in a 
connected region of the address space. 

An important example, where all these issues play a role is the application of 
the Dirac operator onto a source vector in a lexicographic sweep over all lattice 
sites on a node. This is the most time consuming operation during the compu¬ 
tation of propagators. At each site, the source vector, the link matrices and a 
destination vector are needed. In the case of site major they are not contained in 
a connected piece of address space. The pointer to the data of interest therefore 
has to be moved very often and only a single stream is recognized. 

This situation can be improved by copying the necessary fields in lexicographic 
order into temporarily allocated variables at the beginning of the inversion. These 
temporary fields are then used throughout the inversion. Each temporary variable 
now occupies a connected piece of address space which the processor can recognize 
as streams. This philosophy of memory allocation is usually referred to as field 
major (cf. figure 4.2). In order to maintain the user-friendliness of the site 
major approach, the fields are copied back to the site structure ordering after the 
inversion has stopped. This can be done in a negligible amount of time. 
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These ideas have been implemented in the MILC code for the stabilized bi¬ 
conjugate gradient algorithm and the Wilson Dirac operator (dslash_tmp in 
dslash_lean. c and bicgilu_cl in d_bicgilu_cl_cttilde_lean. c). 

The measured performance improvement on a single IBM p690 CPU for dif¬ 
ferent lattice sizes has been summarized in table 4.1. While the achieved speedup 
hrst rises with the lattice volume and then saturates, the overall performance re¬ 
duces with increasing lattice volume. For small lattices, a larger part of the helds 
permanently resides in the cache, which improves performance in general. At 
this point the difference between site major and held major is already sizeable. 
The larger the part of the lattice data residing outside the cache becomes, the 
more the effects of a sensible memory allocation manifests itself. The numbers 
quoted for the speedup are quite impressive and illustrate that it is of utmost 
importance to take into account the target computer architecture at the time of 
the creation of the production code. 

Table 4.1 also contains numbers for the efficiency of the code assuming a peak 
performance of 5.2 GFlop/s for a single IBM processor [?]. The performance 
numbers quoted there are disappointing. On custom designed computers for lat¬ 
tice gauge theory, like the APEMille computer [?,?] for example, efficencies of 
up to 50% are not unusual. When the experiments with the MILC-code on the 
IBM-computers at the HLRN code started, no experience with codes for lattice 
gauge theory, in particular for the application of the Dirac operator with a much 
better performance existed [?]. In the meantime, other groups have achieved bet¬ 
ter performance numbers. For example the NIC/DESY-Zeuthen group [?] have 
measured the performance of their implementation of the SU(3)-Dirac operator 
(Clover improved Wilson action) which uses MPI-based parallelism as well. The 
test was done with a Ib'^-lattice. Although not yet published, they claim to get a 
single-processor performance of 802 MFlop/s on the same computer which is by 
a factor of 2.3 better than the improved version of the MILC code. It would be 
interesting to learn more about their program code and to transfer their ideas to 
the MILC code. 


4.2.4 Miscellaneous other changes 

0(a)-improvement at the boundary 

The 0(a)-improvement term proportional to {ct{go) — 1) (cf. (3.20)) is not im¬ 
plemented in the version of the code that the MILC collaboration shares on the 
web. It is a term that adds to the diagonal of the Dirac matrix at the boundary 
and has been implemented in the routine make_clov. The value of Ct{go) (cf. 
table C.2) has to be specihed in the input hie. 
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L/a 

site major 
performance 
[MFlop/s] 

field major 
performance 
[MFlop/s] 

speedup 

field major 
efficiency 

6 

350 

510 

1.5 

9.8% 

8 

246 

467 

1.9 

9.0% 

12 

187 

349 

1.9 

6.7% 

14 

173 

355 

2.1 

6.8% 

16 

93 

345 

3.7 

6.6% 


Table 4.1: Performance improvement for different lattice sizes after migrating 
from site major to field major, measured on a single IBM p690-CPU with a peak 
performance of 5.2 GFlop/s. 


Large file support for gauge configurations and propagators 

A gauge configuration consists of 

[(L/a)^ X T/a]v x [3 x 3]su(3) x [4]/i x [2]c (4.4) 

real numbers. In the case of double precision arithmetics, a lattice of size {L/afi x 
T/a = 48^ X 96 needs roughly 6.1 GB of memory space. A propagator on the 
other hand consists of 

[{L/a)^ X T/a]v x [3 x 3]su(3) x [4 x 2]Dirac x [2]c (4.5) 

real numbers which, in double precision, for the above lattice size corresponds to 
12.2 GB. Here, the symmetries (3.96) have already been exploited. 

To be able to allocate and store such large address spaces of memory, all 
I/O-routines had to be rewritten to work with a 64Bit address space. 

4.3 Testing the code 

Previous to the production runs, the adopted MILG code had to pass through 
a number of tests which are described in the following. All the tests were done 
on the IBM computer. In particular, the ALPHA-Gollaboration’s lattice gauge 
theory code for the QGD Schrodinger Functional was taken as a reference to test 
all routines of the MILG-code that are involved in the production runs. 

4.3.1 Testing the plaquette 

Up to the pseudo random number generator, the MILG code is based on the same 
gauge update algorithm as the one that is implemented in the ALPHA code. The 
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average plaquette value was taken as a test observable to evaluate the functioning 
of the MILC code gauge-update algorithm. 

One test consisted of determining the normalized average plaquette value 

- 1 ) 

with both, the ALPHA and the MILC code. The results are compiled in the 
following table. 


L^xT 

statistic 

^MILC 
^ P 

■^ALPHA 
^ P 

t-MILC 

Ont 

^ALPHA 

'int 

4^ X 4 
16^ X 32 

1000 

100 

0.62871(37) 

0.59319(3) 

0.62850(32) 

0.59323(3) 

0.051(9) 

0.004(1) 

0.041(6) 

0.003(1) 


In the case of the smaller lattice, the measurements were separated by one heat 
bath and 20 over relaxation steps and in the case of the larger lattice they were 
separated by 100 heat bath steps, each followed by 8 over relaxation steps. The 
integrated auto correlation time Tint is quoted in units of sweeps, not distinguish¬ 
ing between over relaxation and heat bath sweeps. 

As it will be detailed in chapter 6, the production runs with the MILC code 
were done with a geometry of {L/a)^ xT/a = 48^ x 96 at /? = 6.7859. The lattice 
cutoff is already large for this f3 (a/ro ~ 0.0625 —> 6.3GeV)^ and perturbation 
theory for the value of the plaquette should give a rough estimate. 

To evaluate whether the MILC-code produces sensible gauge conhgurations, 
a comparison of the average value of the plaquette to computations in Numeric 
Stochastic Perturbation Theory (NSPT) [?] has been carried out. The result is 
depicted in hgure 4.3. The circles correspond to Monte-Carlo data, generated 
with the MILC code. The corresponding parameters are (100 measurements at 
each value oi (3): 


/5 

6.0 6.1 

6.2 

6.45 

6.7859 

L^xT 

16^ X 32 243 X 40 

243 X 48 

32^ X 64 

48^ X 96 


The solid line (where the error band is barely visible) corresponds to the pla¬ 
quette obtained at 10th order in NSPT. Although the non-perturbative data 
tends to agree with NSPT for larger values of /3, there remains a discrepancy 
at /3 = 6.7859. The MILC code gives [/„ = 0.658732(3), while NSPT yields 
Up = 0.65896(15^ 

Still, the hndings indicate, that the MILC code produces sensible gauge con¬ 
hgurations, since the discrepancy between NSPT and the non-perturbative result 
for the plaquette can be expected to vanish for larger values of (3. 


^The Sommer scale tq = 0.5 fm was employed to convert to physical units (cf. section 5.1.1). 
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Figure 4.3: Average Plaquette U{p) from Numerical Stochastic Petrurbation The¬ 
ory (line with error band) and Monte Carlo simulation (circles). 


4.3.2 Testing the implementation of the correlation func¬ 
tions 


All correlation functions introduced in section 3.6 were implemented in the MILC- 
code and tested against the ALPHA-code. For this test, a PERL-script was writ¬ 
ten (alpha2milc), that converts gauge-conhgurations produced by the ALPHA- 
collaboration’s lattice gauge code into a format, which can be imported by the 
MILC-code. In this way, the correlation functions could be evaluated on the 
same gauge background with both codes. The test was done on a gauge con- 
hguration (16^ x 32-geometry), generated by the ALPHA-code at /3 = 6.0, and 
with seven different values of the hopping parameter with a solver stopping cri¬ 
terion of e = 10“^ (The K-values are the ones that are tabulated in table 5.3 for 
(3 = 6.0). All the correlation functions fA^fp^fv^fr and fp and fy were eval¬ 
uated for the combinations of hopping parameters Ki — Ki, Ki — K 2 , ..., Ki — Ky. 
Averaged over the correlation functions at all times Xq, there was a mean devia¬ 
tion of 0.03%. A maximum deviation of 0.5% was observed for the value of the 
boundary-to-boundary correlation function ky at the largest value of the mass. 

As communicated with the authors of [?,?], where roundoff errors in the same 
correlation functions where studied, the deviations of the magnitude observed 
here lead to systematic errors in the hnal results, e.g. for the decay constant, 
that are much smaller than the expected statistical error. 
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Figure 4.4: The meson mass dependence of the decay constant from a comparative 
test-run between the ALPHA-code (squares) and the MILC-code (circles) aX j3 = 
6.0. AT^eas. = 200. 

4.3.3 Testing the whole setup: A comparative test-run 

This test was carried out in order to simulate the conditions for a production run 
with the MILC code, while having results from the ALPHA code with the same 
simulation parameters as a reference. 

The reference from the ALPHA code was given in terms of 200 measurements 
of the correlation functions necessary to construct the decay constant Fps, i.e. 
/^(xo) and fp. The data was taken from the production runs for [?] at /3 = 6.0. 
As a test observable, the pseudo scalar meson decay constant Fps as a function 
of the inverse pseudo scalar mass 1 /rgrups was studied. The above data contains 
all the correlation functions for seven meson masses. While the light quark mass 
was hxed to the strange quark mass, altogether six heavy quark masses around 
the charm quark mass were simulated for (The hopping parameters are those 
summarized in table 5.3). 

The same number of measurements were then carried out with the MILC 
code, employing double precision arithmetics. The data analysis was done with 
the procedure that will be explained in detail in chapter 6. The dependence of 
the decay constant on the meson mass is plotted in hgure 4.4 which shows, that 
the data agrees within errors. 

Both codes produce compatible results under production conditions. 










Chapter 5 

Simulations in quenched QCD 


The lattice simulations in quenched QCD were done for heavy-light mesons con¬ 
taining a strange quark as the light quark and a quark with a mass around the 
charm quark mass as the heavy quark. The correlation functions for altogether 
six different heavy quark masses were simulated at five different values of the 
lattice spacing in order to allow for a reliable continuum extrapolation. This 
chapter presents the choice of simulation parameters. 


5.1 Simulation parameters 

5.1.1 Setting the scale 

Dimensionful parameters are often given in terms of the Sommer scale tq [?]. 
This is convenient in so far as the attribution of physical units to results in 
the quenched approximation is ambiguous. For example, lattice studies of the 
hadronic spectrum in the quenched theory have revealed inconsistencies in the 
comparison with the experiment of up to 10% [?]. Turning this around, one 
expects different hadronic input to influence the results for observables. 

The Sommer scale has two advantages. On the one hand, [ro/a](/9) has been 
determined precisely over the range of cutoff values [?,?,?] which have been used 
in the scaling study in this thesis. Furthermore, in the absence of dynamical 
quarks, tq is only affected by cutoff effects of order [?]. 

When using physical input from experiments, the conversion will be done by 
setting To = 0.5 fm. For this value of the reference scale, the decay constant 
of the iF-meson takes its physical value tqFk = 0.405(5) in quenched lattice 
simulations [?]. One therefore often refers to this choice as setting the scale with 
the Kaon decay constant. Setting the scale with the nucleon mass instead, the 
reference scale would have to be Tq = 0.55 fm to match the experimental value. 
This corresponds to a 10% shift in the scale. An estimate for the size of the 
ambiguity for the observables determined in this work will be given in chapter 
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L/a 

T/a 

L/tq 

a[fm 

^meas 

^Update 


Pi 

= 6.0 

16 

32 

2.98 

0.093 

380 

100 

8 

P2 

= 6.1 

24 

40 

3.79 

0.079 

201 

100 

12 

Ps 

= 6.2 

24 

48 

3.26 

0.068 

251 

100 

12 

Pa 

= 6.45 

32 

64 

3.06 

0.048 

289 

100 

12 

Pb 

= 6.7859 

48 

96 

3.00 

0.031 

150 

50 

24 


Table 5.1: Lattice geometries and simulation parameters used in the simulations 
for the scaling study. 


6.2.4. 

A parameterization of the ratio a/r^ has been obtained as a function of the 
bare coupling /3 = 6 /( 7 q in [?], 



exp {-1.6804 - 1.7331(/3 - 6) + 0.7849(/? - 6)^ - 0.4428(/? - 6)3} , 


for 5.7 < /3 < 6.92. 


(5,1) 

The accuracy of this formula ranges from 0.5% at low values of jS to 1% at large 
values of jS. 


5.1.2 Parameters for the scaling study 

For the scaling study, simulations were done for hve different lattices of ap¬ 
proximately constant physical size L/tq 3, but decreasing lattice spacing. 
The values /3i ,..., Pi were taken over from previous work within the ALPHA- 
collaboration [?]. The value of was obtained with (5.1). The corresponding 
lattice geometries are given in table 5.1. 

As the time-extent of the lattice, T k, 2L was chosen. The experience from 
earlier simulations with the Schrodinger Functional, e.g. in [?, ?] showed, that 
with this asymmetric geometry, the plateaus of the effective masses (3.76) and 
decay constants (3.78) are well pronounced. A numerical check for pseudo scalar 
mesons in [?] revealed, that with a lattice extent oi L/tq k, 3 fm, T = 2L, 
hnite volume effects can be neglected when simulating for mesons containing a 
strange quark as the light quark. The corresponding vector mesons are heavier 
and therefore have a shorter Compton wavelength. It is concluded, that finite 
volume effects for this channel are also negligible. 

A peculiar choice for the lattice geometry had to be made at P 2 = 6.1, 
where the hardware and memory conhguration of the APEMille computers at 
NIC/DESY-Zeuthen [?] did not allow to follow the relation T = 2L. 
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5.1.3 Improvement and renormalization constants 


Most of the improvement and renormalization constants needed for the derivation 
of simulation parameters and for the data analysis have already been introduced 
in chapter 3. Where available, only non-perturbatively determined improvement 
and renormalization constants were used. The parameterizations of all constants 
in terms of the bare coupling qq or f3 = 6 /Qq are given in appendix C, in the 
tables C.l and C.2. 

Two additional renormalization constants, namely Z [?] and Zm [?] will be 
needed. On the one hand. 


rrii 

rriq^i 


(5.2) 


relates the bare current quark mass (PCAC) of some flavor i to the subtracted 
bare quark mass while 


Zm{9o) 


Mi Za^Qo) 
mi(/r) Zp{go,fi) 


(5.3) 


on the other hand that relates the bare current quark mass to the renormalization 
group invariant mass Mj. 

In [?], the non-perturbative parameterization of Za{P) [?] together with sim¬ 
ulations for the ratio Mi/rfii{fi = (2Linax)~^) = 1.157(12) and for Zp{go,ii = 
( 2 Linax)~^) at five values of /3 were combined to give a parameterization for 
ZM^go) in the range 6.0 < P = 6/g^ < 6.5. The scale dependence of the pa¬ 
rameters was computed in the mass independent Schrodinger Functional scheme, 
where the renormalization scale is given by the hnite box size, L = l//i. In this 
particular case, /i = (2L max )~^ ~ (1.436ro)“^ was used. 

With further data from [?], it is possible to derive a parameterization of 
ZM{go) for a larger range of jS which will be needed in the continuum extrapolation 
for the renormalization group invariant charm quark mass (cf. table 5.1). 

First, one has to find an extended parameterization of Zp{go, (2Lmax)~^)- This 
factor has been determined for several values of /?, once along lines of constant 
physics defined by the fixed volume Umax and once along an equivalent (in the 
continuum) condition defined by the renormalized coupling in the Schrodinger 
Functional [?], Zp{gQ, ( 2 L)“^)|g 2 (i/i )=343 [?]. The corresponding data is summa¬ 
rized in table 5.2. The two definitions of lines of constant physics differ by 0{a^), 
since the simulations have been done in the improved theory. 

The two data sets are compatible, as can be seen in hgure 5.1. The diamonds 
represent the data along constant volume and the circles along constant renor¬ 
malized coupling. The parameterizations for Zp{gQ, (2Lniax)~^) [?] and a new 
parameterization for the data of Zp{gQ, ( 2 L)“^)|g 2 ( 4 /j ^)=3 43 have been added to 
the plot (dash-dotted and dashed line). The latter parameterization describes 
both the data sets very well in the range 6.0 < j3 < 7.0 and the combination 
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hxed 2Lmax ~ 

1.436ro 




6.0 6.0219 

6.1628 

6.2885 

6.4956 

Zpi^gOy (2Tniax) ) 

0.5253(26) 0.5218(16) 

0.5177(19) 

0.5179(23) 

0.5157(19) 


hxed g^{l/L) = 

3.480(13) 




6.257 6.476 

6.799 

7.026 


Zp{go, (2L)-1) 

0.5179(19) 0.5143(23) 

0.5133(19) 

0.5137(27) 



Table 5.2: Non-perturbative data for Zp{gQ, (2Lmax) ^)- 


of both data sets into one parameterization is thus justihed. For a polynomial 
ansatz one obtains 

Zp{j3, (2Lj^ax)”^) = 0.5228 - 0.0231(/3 - 6) + 0.0142(/3 - 6)^. (5.4) 

It describes the data with a maximal deviation of 0.5% and is plotted in hgure 
5.1 as the solid line. 

In a second step, the data for Zp{gQ, (2Lmax)~^) can be combined with the 
ratio Mj/mi((2Lrnax)~^) = 1.157(12) and with the non-perturbatively determined 
values of Za{(3) to give ZM^go) which then can be parameterized as 

Zm{P) = 1.754 + 0.27(/3 - 6) - 0.10(/3 - 6)1 (5.5) 

This parameterization describes the data within 0.5% accuracy. The error of 
ZA{go) and of the ratio Mj/mj((2Lmax)~^) was taken into account in the error 
analysis. 

5.1.4 Hopping parameters and stopping criteria 

In order to obtain correlation functions for mesons containing a strange quark 
and for various heavy quark masses in the range of charm, propagators for seven 
different hopping parameters were determined in the Monte-Carlo simulation. 
They were chosen such, that the lightest quark mass corresponded to the strange 
quark mass. A second hopping parameter was tuned to simulate at the charm 
quark mass and the remaining hve parameters were distributed homogeneously 
over the mass range between the strange quark mass and the lattice cut-off. 


Kcrit - the critical hopping parameter 

Kcrit has been determined in [?] in quenched QCD for a large range of /3-values 
(c.f. table 1 therein). For /3 = 6.0 and 6.2, the values could be taken over 
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Figure 5.1: The renormalization factor Zpjgn, l/2L m »^)- The diamonds cor¬ 
respond to non-perturbative data at hxed 2Linax = 1.436ro and the circles 
correspond to non-perturbative data obtained from the simulations at hxed 
f{l/L) = 3.48. 
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Figure 5.2: Quadratic interpolation for HcritiP)- Only the data at the hlled squares 
was included into the £t. 


directly. In collaboration with the authors of [?], the critical hopping parameters 
for (3 = 6.1, 6.45, and 6.7859 were determined from quadratic interpolations in /?, 
of which the case (3 = 6.7859 has been illustrated representatively in hgure 5.2. 


Kg - the hopping parameter corresponding to the strange quark 


The hopping parameters for the strange quark for (3i — (3^ were obtained, using 
previous work by the ALPHA- and UKQCD-collaboration [?]. There, the sum 
of the renormalization group invariant strange and light quark mass Mg + M 
with M = \{Mu + Md) has been determined in a quenched simulation with 
0(a)-improved Wilson fermions. The basic idea in this work was to exploit the 
PCAC-relation 


Mg + M = Zu^ml (5.6) 

between the renormalization group invariant quark masses of the strange quark. 
Mg, the average light quark mass M = |(M„ -f M^) and the mass of the corre¬ 
sponding meson, mx. Zyi, as detailed in the previous section, relates the bare 
current quark mass to the renormalization group invariant quark mass. Fk is the 
Kaon decay constant and Gk denotes the vacuum-to-JL matrix element of the 
pseudo scalar density, which was determined in that work [?]. 
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(a/ro)2 

Figure 5.3: Linear extrapolation in (a/ro)^ for with la-error-band. Only 

the data at the hlled squares was included into the ht. The circle represents the 
extrapolated value at /? = 6.7859. 


The computation used rorriK = 1.5736 [?] as experimental input. The follow¬ 
ing table summarizes the results for aX (3 = 6,6.1,6.2,6.45 obtained with 

0(a)-improved Wilson fermions. The value at /? = 6.7859 has been extrapolated 
linearly in {a/ro)"^ based on the numerical data at /? = 6.1, 6.2, 6.45, as is depicted 
in hgure (5.3). 


/9 

6 

6.1 

6.2 

6.45 

6.7859 

ZmFk 

GKro 

0.1939(30) 

0.2077(28) 

0.2160(30) 

0.2205(46) 

0.2268(57) 


Mg can be extracted from (5.6) by using the ratio Mg/M = 24.4 ± 1.5 from chiral 
perturbation theory [?]. Together with the dehnition of the subtracted bare quark 
mass (3.46), this dehnes the hopping parameter Kg of the strange quark as the 
solution of the quadratic equation in atUg^g = 1 — «:“(^(5fo)), 


ro{Mg + M)=roMg{l + —) 



ZM^ClTflqgiX -|- bmOiTnqg^. 


(5.7) 
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(a/ro)2 

Figure 5.4: Linear extrapolation in [a/r^Y for ^o^c- Only the filled squares 
entered the £t. The circle represents the extrapolated value at /3 = 6.7859. 

Kc - the hopping parameter corresponding to the charm quark 

The hopping parameters at /3i... /34 for the charm quark have been determined 
in the computation of the renormalization group invariant mass [?]. With the 
values for Me at finite lattice spacing (/3i — (3Yj listed in the following table, one 
can extrapolate linearly in (a/ro)^ to obtain a first guess at (3^. The extrapolation 
is illustrated in (c.f. figure 5.4). 


p 6 6.1 

6.2 

6.45 

6.7859 

roMe 3.224(41) 3.479(43) 

3.711(47) 

3.975(53) 

4.277(85) 


As in the case of the strange quark hopping parameter, Kc can again be obtained 
by solving the quadratic equation (5.7). 

Five additional hopping parameters 

Five additional hopping parameters were guessed such that a roughly uniform 
distribution of quark masses between the strange quark and half the 5-quark 
mass for /3i,... ,/34 and 4.5 GeV in the case {3^ was achieved. In particular, the 
parameters were determined from a linear inter- resp. extrapolation with respect 
to the pseudo scalar masses romx = 1.5735 and romo^ = 4.988 [?]. 
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6 

6.1 

6.2 

6.45 

6.7859 

^crit 

0.135196 

0.135496 

0.135795 

0.135701 

0.135120 

Ki 

0.134108 

0.134548 

0.134959 

0.135124 

0.134739 


0.128790 

0.130750 

0.131510 

0.132690 

0.132440 

^3 

0.123010 

0.125870 

0.127470 

0.130030 

0.130253 

Ki 

0.119053 

0.122490 

0.124637 

0.128131 

0.128439 

K5 

0.115440 

0.119370 

0.122000 

0.126330 

0.126774 

Kq 

0.112320 

0.116640 

0.119680 

0.124730 

0.123571 

Kf 

0.109270 

0.113960 

0.117370 

0.123120 

0.117625 


Table 5.3: Summary of all hopping parameters. The shaded helds indicate the 
value of the hopping parameters corresponding approximately to the strange 
quark and to the charm quark. 


The solver stopping criterion e 

Heavy-light correlation functions decay exponentially with the meson mass in 
time. Thus, large components contribute in the heavy quark propagator for 
initial times and small components for large times Xq. If one takes for example the 
heavy-light correlation function /^(a^o)) that decays by 28 orders of magnitude 
between xq = a and Xq = T — a for the cases with a very large heavy quark 
mass ((L/a)^ x (T/a) = 48^ x 96, (3 = 6.7859, and kj as in table 5.3 ). In 
its computation, sums over numbers of different order of magnitude have to be 
evaluated, thereby possibly introducing roundoff errors if the number precision is 
not sufficient. 

Adjusting the stopping criterion of the inverter is a crucial task in the prepa¬ 
ration of the production runs. A stricter value for e causes the inverter to iterate 
longer. But the stopping criterion is also connected with the precision to which 
the observables constructed from quark propagators are determined. The runs 
at /9 = 6.0, 6.1, 6.2 and 6.45 were all done in single precision and the value of 
e = 10“^ was taken for all hopping parameters. The run with the MILC code at 
(3 = 6.7859 used double precision arithmetics. In order to hnd out the optimal 
stopping criteria for this production run, the quark propagators for all seven val¬ 
ues of the hopping parameter were computed on a gauge background, once with 
e = 10“^ and once with e = 10“^®. All meson correlation functions dehned in 
chapter 3.4, were then computed for all combinations of the hopping parameters 
Ki — K 2 , Ki — ks, ..., Ki — Kj. Then, the relative deviations of correlation functions 
were computed. It turned out, that especially the light quark propagators are 
not sensible to the changed stopping criterion (cf. figure 5.5). The maximum 
deviation for the combinations of hopping parameters — ^ 2 ; • • •, ~ ^4 was 

less than 0.01% and the mean deviation of all correlation functions at all values 
of Xq was better than 0.002%. In contrast, for the combinations ki — ks, ki — kq 
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Figure 5.5: Upper plot: Relative deviation between correlation functions evalu¬ 
ated once with a stopping criterion of e = 10“^ and once with e = 10“^®. Lower 
plot: Relative deviation between e = 10“^^ and e = 10“^^ (plus signs) and be¬ 
tween e = 10“^^ and e = (circles). Each plot shows the relative deviation for 

the correlation functions fA^xo), fp{xo), kv{xo) and kpixo) for xq/q = 1... 95. 
Missing markers indicate vanishing relative deviation. 
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and Ki — Kj the correlation functions for the two different stopping criteria dif¬ 
fered by orders of magnitude for times Xq/q > 50. In a second test, the meson 
correlation functions for the combination of hopping parameters ki — were 
computed with stopping criteria e = 10“^^, 10“^"^ and 10“^^. The lower plot in 
figure 5.5 shows the result. While the correlation functions deviate sizably when 
changing the solver precision from e = 10 “^^ to e = 10 “^^ (crosses), the change is 
small between the two solver precisions e = and e = 10“^^ (circles). Thus, 

a conservative solver residual of e = 10 “^® for the large quark mass associated to 
^5 was used. The same stopping criterion was applied to kq and ^ 7 . However, 
roundoff errors could not be ruled out in these cases. 


5.2 The simulations 

The simulations at /3 = 6.0, 6.1, 6.2 and 6.45 were all done on the APEMille 
computer at NIC/DESY-Zeuthen [?] with the ALPHA-collaboration’s TAO-code, 
while the simulation at /? = 6.7859 was done with the MILC-code at the HLRN 
[?]. The data of the former runs were already available when this project started 
and could be taken over. 

Beginning with a cold gauge configuration, 0(100) update sweeps were done 
at each value of jS to guarantee thermalization. The average plaquette value 
(4.6) was also monitored but took a stable mean earlier at every (5. The runs 
on the APEMille computer were not so time consuming and a large number of 
intermediate updates, followed by 0{L/2a) over relaxation steps was done at all 
values of j3. The number of intermediate updates was reduced for the run at 
(3 = 6.7859, but no correlation effects in the data were observed. 

After the thermalization phase, the conhgurations for the production were 
generated. 

The runs on the IBM computers at the HLRN 

The generation of gauge conhgurations was carried out with a parallelization over 
64 CPUs where the 50 update steps take 3.4 hours. 

To save CPU-time, the generation of the held conhgurations was done using 
the code with single precision arithmetics.The gauge conhgurations obtained in 
this way were then converted to double precision arithmetics using a PERL-script. 

During the computation of the correlation functions with double precision 
arithmetics, a total of 114 GB of memory had to be allocated: 6 GB for the gauge 
held, 56 GB for the storage of one color and two spin components of the seven 
quark propagators and the rest for auxiliary helds, necessary for the biconjugate 
gradient routine. This large amount of memory assigned to the auxiliary helds 
was due to the re-shuffling which was necessary for changing from site-major to 
field-major (cf. 4.2.3). 
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As the IBM servers at the HLRN are divided into units of 8 CPUs (LPAR) 
with a shared memory of at least 64 GByte, it was expected, that the code per¬ 
forms properly with a parallelization over 64 CPUs. In this case, the program 
allocates 57GB of RAM on each LPAR. However, the code then ran very instable. 
Fluctuations of the total runtime by a factor of two were observed. The project 
consultants at the HLRN found out, that memory intensive system software (e.g. 
the operating system itself and software connected to the CPFS hie system man¬ 
agement) consumed a considerable part of the shared memory. This caused the 
production code to swap, i.e., parts of the allocated memory temporarily had to 
be stored to disk in order to free dynamic memory, thereby slowing down the pro¬ 
gram considerably. In order to avoid these losings in performance, the production 
was reorganized to use 128 CPUs instead. This unfortunately prolonged waiting 
times in the Queue, because larger numbers of CPUs were free less frequently. 

With a parallelization with 128 CPUs, the computation of the propagators 
for the seven hopping parameters given in table 5.3 took 9.5 hours. Loading the 
gauge conhguration in the beginning, contracting the propagators to compute the 
correlation functions and saving propagators took a negligible amount of time. 

Application for CPU-time 

In the beginning, the plan was to produce a statistic of 200 measurements of all 
the forward and backward correlation functions at (5^. 

The thermalization and most of the initial tests could be done at the HLRN 
before the official accounting started and was for free. For the production runs 
however, applications for CPU-time had to be written. 

At the HLRN, the CPU-time is measured in NPL [Norddeutsche Parallel- 
rechner-Leistungseinheit). 1 NPL is dehned as 1 hour wall-clock time on 32 IBM 
p690 CPU’s. 

At the time of writing of the hrst application, some of the improvements of 
the code were not yet implemented and one point in the Monte-Carlo-History was 
estimated to cost 128 NPL. On this basis, an application for 25600 NPL for one 
year was written, of which 10000 NPL were granted for the period 04-2003 until 
03-2004. Another 10000 NPL were granted after an application for an upgrade 
in 09-2003. 

The progress during the hrst six months of the production runs for (3 = 
6.7859 was very slow because the computing facilities did not provide the expected 
services. On the one hand, the computer at the HLRN had a lot of down times 
and the code did not perform very well, as detailed above. Furthermore, a fair 
queueing system, that prefers massively parallel applications from jobs with a 
small number of CPUs was introduced only after the hrst six months of the 
production. 

Most of the problems were hnally resolved after nine months when the project 
was granted an extension of another three months. At this point, the generation of 
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one gauge configuration cost 7.8 NPL and the measurement of the corresponding 
forward and backward correlation functions roughly 80 NPL. 

The project hnished with a statistic of 151 measurements at /3 = 6.7859. The 
forward propagators and the backward propagators for the hopping parameters 
Ki and K 3 , which approximately correspond to the mass of the strange quark and 
the charm quark respectively, were stored on tape for 80 measurements. They 
can be restored for further usage. 



Chapter 6 
Data analysis 


This chapter addresses the analysis of the data which were obtained from the 
Monte-Carlo simulations of lattice QCD in the quenched approximation. The 
continuum limit was taken for the decay constant and the mass splitting of the 
Ds*^-meson and the renormalization group invariant charm quark mass. The 
combined analysis of the simulation results for a number of heavy quark masses 
around the charm quark mass, together with predictions and simulation results 
from HQET, allowed to determine the decay constant of the Bs*^-meson and 
the corresponding mass splitting at the physical point from an interpolation in 
the meson mass. Furthermore, estimates for the magnitude of the first spin- and 
flavor-symmetry breaking terms in the l/Mg-expansion in HQET could be given. 


6.1 Data analysis - general remarks 

The results and errors from the simulations at all values of (3 have been obtained 
from statistical samples of the primary quantities 


fo{xo),go(T - xq) for O = A, V, P, T and 


for 0 = P,V, 


( 6 . 1 ) 


over the whole range of Xq G [a,T — a], using the jackknife method. The correla¬ 
tion functions for the forward and the backward direction /o(a^o) and go(T — xq) 
have been averaged at each step of the Monte-Carlo history^, 

fo{xo) = ^{fo{xo) + go{T - xo)). (6.2) 

^This is possible only for vanishing background field in the Schrodinger Functional. 


63 



64 


CHAPTER 6. DATA ANALYSIS 


Using the improvement and renormalization constants introduced in section 5.1.3, 
the secondary quantities 

Fps, Fv, Fps/Fv, mps, my, mpg - my and 

(6.3) 

Yps Yy R Am 
Cps ’ Cv ’ C'ps/v’ C'spin 

were constructed from them following the definitions in section 3.6.2 and section 
2.3. No autocorrelation in the data was observed. 

The renormalization and improvement coefficients are only known up to sta¬ 
tistical and systematic errors from their determination in lattice simulations. In 
order to take this error properly into account, a statistical sample of all constants 
with a Gaussian distribution around their mean and with the width defined by 
their error was generated, using a pseudo random number generator (randn in 
MATLAB). These samples were handed over to the jackknife routine and then 
treated in the same way as the Monte-Carlo data of the primary observables. 

The error due to perturbation theory in the conversion functions Cy^^Mq/ 
as listed in table 2.3 has been taken into account by error propagation. 

6.1.1 Plateaus 

First, jackknife results were determined for the effective mass 

"•2,(^0+ i) = (6,4) 

Figure 6.1 shows the results for the combination ki — of hopping parameters 
at /ds = 6.7859. As expected from (3.76), the effective mass exhibits a plateau for 
intermediate times, and contributions from excited states of mass A and glueballs 
of mass me for small and large times xq, respectively. 

The meson mass my can be extracted as the average over the plateau. In order 
to keep the systematic errors in my due to contaminations by excited states under 
control, the time interval, where their relative contributions are below a chosen 
threshold was determined. The thresholds that have been used are given in table 
6.1. Thus, the statistical error will always exceeds the systematic error. The time 
interval, or plateau range, can be determined by means of the following iterative 
procedure: 

One first subjectively chooses a sensible plateau range and subtracts the av¬ 
erage over the plateau from the data. 

1.) From the logarithm of this data, estimates for the contributions (cf. (3.76)) 

2 sinh(aA/2)?7®e“'^°^, 2 sinh(amG/2)?7xe“*^^“*°)™° (6.5) 

to the effective mass can be obtained in terms of linear fits to the time 
dependence for small and large times. The fit ranges have to be chosen 
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xo/a xo/a 


Figure 6.1: Plot of the effective mass and the decay constant Fx(xo) at 
(3 = 6.7859 for the combination of hopping parameters ki and 


romps 

romv 

’^oFps 

roFy 

0.5% 

0.7% 

0.5% 

0.7% 


Table 6.1: Thresholds for the accepted contribution of excited states to the 
plateau range. 
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subjectively. This procedure is illustrated in figure 6.2, again for ki — 
% at /3 = 6.7859. The data sometimes does not exhibit a clear linear 
behavior and the ht range cannot be chosen without ambiguities. Therefore, 
the extracted glueball masses and mass gaps can only be interpreted as 
estimates, which however suffices for the purposes here. 

2. ) The sum of the relative contributions of excited meson states and glueballs 

to the plateau as a function of the time xq is shown in hgure 6.3. The 
new plateau range is dehned as the time interval, where this contribution 
is below the threshold given in table 6 . 1 . 

3. ) The procedure can be repeated with the newly dehned plateau average, 

until the plateau range is stable, which usually occurs after one or two 
iterations. 

One hnally obtains the meson mass mx as the average over all the plateau 
ranges listed in table D.l in the appendix. They typically extend over Irg to 1.5ro 
and their position approximately scales with j3. The estimates for the glueball 
mass have been summarized in table 6.2. The values given there for each value 
of jS are the averages over the hts obtained from the six hopping parameter 
combinations ki — K 2 ... ki — K 7 . The estimates for the gap energy are tabulated 
in table D.2. 

With the plateau ranges of mx at hand, one now has to repeat the procedure 
for the decay constant 

Fx(xo) = - 2 Zo(mxL 3 )-V 2 ePo-T/ 2 )mx^_ 

Here, the contributions of the excited states to the plateau are of the form 

Representative plots of the procedure are shown next to the ones for the effective 
mass in the hgures 6.1 to 6.3. 

Although feasible in principle, a continuum limit for the glueball mass has not 
been taken because of the uncertainties in the determination of the corresponding 
ht ranges. But even at hnite lattice spacing, the data is mostly compatible 
with the glueball mass of the 0 ’''’''-glueball as obtained in the dedicated lattice 
simulations and summarized in table 6.3. 

For the determination of the charm quark mass using the PCAC-relation, no 
plateau average over the involved correlation functions has been taken, as the 
statistical errors turned out to be very small. Instead, the value at Xq = 2/3T 
was used, where the contamination by excited states was reasonably small. 


(6.7) 
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Figure 6.2: Fits to the logarithm of the subtracted effective mass and the sub¬ 
tracted decay constant. The slope of the linear hts give an estimate for the mass 
gap A (dashed line) and the glueball mass niQ (dash-dotted line). The dotted 
lines indicate the corresponding £t range in each case. 



r 0"lG,Fps 

romG,Fv 

^O^G.mfll 


A 

5.9 

- 

5 

4 


5.3 

4.9 

5 

4 


4.5 

4.5 

4.1 

3.7 


4.7 

4.4 

4.0 

3.3 

/55 

4.38 

3.4 

3.5 

3.5 


Table 6.2: Estimates for the glueball mass. In the case r’omG,Fv A; the data 
was too noisy to allow for a sensible ht. 
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xo/a xo/a 



xo/a xo/a 

Figure 6.3: Estimated systematic contribution to the plateau due to excited 
meson states (dashed line), due to glueballs (dash-dotted line) and the sum of 
both contributions (solid line). The plateau average is taken over the range, 
where both contributions are below the threshold (indicated by the bold dashed 
line). 


Collab. 


romo++ 

Morningstar et. al. 

[?] 

4.21(11)(4) 

Bah et. al. 

[?] 

4.33(10) 

Teper 

[?] 

4.35(11) 

UKQCD 

[?] 

4.05(16) 

Niedermayer et. al. 

[?] 

4.12(21) 


Table 6.3: Lattice data for the lowest glueball mass (0'*’+) in the continuum with 
(combined) statistical and systematic errors. 
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6.1.2 Interpolation in the meson mass 

The observables at different quark masses but at constant jS have been evaluated 
on the same set of gauge configurations and are therefore correlated. Indeed, 
the corresponding normalized covariance matrix of the observables which was 
determined in the jackknife procedure has large off-diagonal elements. It was 
checked, that this does not affect the error estimation in the data interpolation. 

Interpolation to lines of constant physics 

The pseudo scalar and vector meson masses computed at each lattice spacing 
do not coincide exactly (cf. figure 6.4 and table D.3 in the appendix). Thus, 
for a continuum limit along a line of constant physics, the secondary quantities 
(6.3) were interpolated to common values. For the decay constant and the mass 
splitting, this was done linearly in the inverse meson mass mps or my, whereas 
in the case of the renormalization group invariant charm quark mass, linearly in 
the meson mass. These parameterizations are not only suggested by HQET, but 
are also supported by visual examination of the mass dependence of all secondary 
quantities. 

Figure 6.4 shows the choice of meson masses given in table 6.4 for the pseudo 
scalar and the vector meson channel, to which the secondary observables were 
interpolated. Most of the masses were chosen such as to lie in the vicinity of the 
simulated ones in order to reduce the error introduced by the interpolation. In 
one case, the interpolation was performed to the experimentally known masses [?] 
of the Dg and D* meson respectively. 

The data obtained for kq and K 7 at (3^ = 6.7859 were discarded for the whole 
analysis, since roundoff effects were observed for the heavy quark masses simu¬ 
lated for, which could be ruled out reliably only for (cf. 5.1.4). 

Interpolation to the 5-quark 

Prior to the interpolation between the region of the charm quark mass and the 
static limit, the continuum extrapolation for the desired observables at the masses 
given in table 6.4 was carried out. This ordering of the procedure is suggested by 
the expectation, that discretization errors are quark mass dependent. Since all 
the observables have been determined in the 0 (a)-improved theory, one would 
naively assume, that they approach the continuum limit approximately linearly in 
a?. In [?], deviations from this scaling behavior for heavy quarks were observed. 
A comparison of lattice results in the finite volume 0(a)-improved Schrodinger 
Functional at non-vanishing lattice spacing with results from perturbation theory 
at zero lattice spacing revealed, that the expected scaling behavior breaks down 
for heavy quark with (arug®)^ > 0.2. Thus, the bound 


aMq < 0.64. 


(6.8) 
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i 

1 

2 

3 

4 

5 

6 

rom|,s 

3.768 

4.327 

4.987 

5.653 

6.211 

6.560 

aMq < 0.64 

^2 — Pb 

h — Pb 

/ds — Pb 

/ds — Pb 

/ds — Pb 

/ds — Pb 

romV 

4.210 

4.660 

5.363 

5.920 

6.280 

6.550 

aMq < 0.64 

^2 — Pb 

P2 — Pb 

/ds — Pb 

/ds — Pb 

/ds — Pb 

/ds — Pb 


Table 6.4: Meson masses to which all observables were interpolated before taking 
the continuum limit. The masses in the third column [i = 3) correspond to the 
experimentally known values of the Ds and D* meson respectively [?]. In addition, 
the lattices which will be included into the continuum extrapolation are given in 
terms of the corresponding coupling constant (cf. section 6.1.2). 
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Figure 6.4: Simulated inverse pseudo scalar masses for the combinations of hop¬ 
ping parameters ki — K 2 , ■■■, Hi — (circles). The dashed lines indicate the 
inverse masses to which the observables have been interpolated at each value of 
/3 before the continuum limit has been taken. 
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has to be fulfilled by all the data that enter the analysis. Table 6.4 shows the 
meson masses together with the range of /5-values that are compatible with this 
bound. 

The interpolations were carried in the meson mass. This allowed to localize the 
physical point of the Bg-meson exactly, because its mass is known very precisely 
from experiment (mB^ = 5.36966(24) GeV [?]). In units of the Sommer-scale, the 
mass translates to rouiB^ = 13.6056(6). 

In HQET, one expects that the heavy quark mass differs from the heavy-light 
meson mass only by the binding energy A and terms of 0{1/Mq) (cf. (2.43)). 
Thus, the functional dependence of the observables on the meson mass is com¬ 
patible with the fit-ansatz 


Ul 02 

do H-1-2— 

romps rorups 


+ ... . 


(6,9) 


The coefficients a* can be taken as estimates for the magnitude of the contribu¬ 
tions to heavy-light observables from higher orders in the HQET expansion. 


6.2 The Dg- and the D*-meson and the c-quark 
mass 

In this section, a precise determination of the observables 

Ed,, Ed*, Ed./Fd,* andro(mDj - moj (6.10) 


is presented. 

Simulations with a continuum extrapolation for the decay constant from 
rather coarse lattice resolutions have been done in [?,?,?,?] and more recent 
studies with an extended approach to the continuum limit, 0(a)-improved Wil¬ 
son fermions and also partly with dynamical quarks have been carried out in 

A recent summary of lattice data for the Dg-system can be 
found for example in [?]. Some representative results are summarized in table 
6.5. 

6.2.1 The decay constants Fd^ and Fd* and the ratio Fd^/Fd* 

After interpolating the data for the decay constants and Fo* at each finite 
lattice spacing (cf. table D.3) to the physical point, given in terms of the mass 
^omDs = 4.987 [?], one obtains the corresponding values in the continuum from 
a linear extrapolation in (a/ro)^. The extrapolations for and Fb^/Fdj are 
shown in figure 6.5. For Fd„, the fit has y^/d.o.f. = 1.9, for F^j y^/d.o.f. = 1.6 
and for Fd^^/Fd* a value of x^/d.o.f. = 0.7. The final results in the continuum 


are 
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Figure 6.5: Continuum limit for the decay constant Fd^ and the ratio Fd^/Fd*. 

roFc, = 0.572(18) , 
toFd* = 0.605(47) , 

Fd./Fdj = 0.939(61) . 

Using ro = 0.5 fm, the results for the decay constants translate to 

Fds = 226(7)MeV and Fd* = 239(18)MeV. The error is the combination of 
statistical and systematic contributions within the quenched approximation. 

Following the arguments of section 6.1.2, the data from the simulations at 
Pi = 6.0 and P 2 = 6.1 has not been included into the hts. However, the stability 
of the extrapolation has been tested by studying the change in observables under 
the inclusion of P 2 = 6.1. This resulted in roFD,,=0.585(15) with y^/d.o.f. = 1.9, 
^"oFd*= 0.658(38) with = 2.4 and Fd,,/Fd*= 0.885(49) with y^/d.o.f. = 1.5. 
Thus, the observables change by 2%, 8% and 6% respectively, and the results 
from both fit ranges agree within errors. 

The data at P 2 — Pa has previously been used to determine the pseudo scalar 
decay constant in the continuum with a hnal value of Fd^ = 252(9)MeV [?], 
where the scale was also set with the Kaon decay constant. This value differs 
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Reference 


F^f "/MeV N, 

FD./MeV 

F^r^/MeV 

scale 

setting 

ALPHA 

[?] 

252(9) 




ro 

Becirevic et. al. 

[?] 

231(12)(«) 



272(16)(«) 

rrip 

Bowler et. al. 

[?] 

229(3) (/g) 



264(10)(/‘a 

U 

CP-PACS 

[?] 

260(l)(/g) 

2 

267(13)(«;) 


rrip 

de Divitiis 

[?] 

240(5)(5) 




ro 

MILC 

[?] 

223(5)(i!|) 

2 

241(6)(/») 


U 

UKQCD 

[?] 

229(3) («") 




U 

Wingate et. al. 

[?] 


3 

290(7)(41) 


T 

Ryan {world av.) 

[?] 

230(15) 


250(30) 




Table 6.5: Results for the decay constant and Fd* from other groups. The 
errors are statistical and systematic. 


from the one obtained here by 26 MeV and the two results are not compatible 
within errors. 

Table 6.5 summarizes some recent values for the decay constants from quenched 
as well as from dynamical simulations by other groups. There are indications, 
that the effect of unquenching is a shift in the decay constant to about 10% higher 
values [?,?]. However, results with dynamical quarks still suffer from large sta¬ 
tistical and systematic uncertainties so that a reliable estimate of the quenching 
error is not yet possible. 

Only the decay constant for the Dg-meson has been determined in an ex¬ 
periment. Two recent measurements quote F“^ = 285(19)(40)MeV [?] and 
FqP = 280(19)(28)(34)MeV [?], where the hrst error is statistical and the sec¬ 
ond systematic. In the latter case, the third error is also systematic due to an 
uncertainty in branching ratios. The present world average from experimental 
determinations is F^^'’ = 267(33)MeV [?]. 

Also QCD sum rules have been used to determine Fd^. In [?] a value of 
Fd„ = 235 (24)MeV has been suggested. 

6.2.2 The mass splitting ro(mD* — 

In the same way as for the decay constant, after interpolating the data to the 
physical point given by the mass of the Dg-meson, the continuum extrapolation 
can be carried out for the mass splitting ro(mD* — (cf. hgure 6.6). The 

resulting value is 


ro(mD* - moj 


0.345(23), 
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(a/ro)^ 

Figure 6.6: Continuum limit for the mass splitting. The diamond represents the 
experimental value. Only the data at the hlled squares entered the £t. 

The linear fit has x^/d-O-f- = 0.1. Converted to physical units with tq = 0.5 fm, 
the value is (rriD* — moj = 136.0(92) MeV. The result is compatible with the 
current experimental data taken from [?], (uid* — = 143.8(4)MeV. 

In a previous lattice study, where the same techniques were applied but the 
continuum extrapolation was made from coarser lattices [?], a value of (ttid* — 
tudJ = 122(14)MeV was quoted. Another computation [?] gave (mr,* — tudJ = 
97(12)MeV as the final result. This result was obtained on rather coarse lattices 
and the discrepancy with respect to experiment was ascribed to cutoff effects. 


6.2.3 The renormalization group invariant charm quark 
mass Me 

The renormalization factor Zm(5'o)) which has been introdneed in section 5.1.3, 
can be used to relate the two definitions of renormalized quark masses (3.47) and 
(3.81) to the renormalization group invariant quark masses. 

On the one hand, employing the dehnition of the PCAC-mass (3.80), the 
charm qnark mass can be extracted from the correlation functions containing 
non-degenerate quarks. 


'^O-^ejmsc 


ZM\2romsc [l + (^a - bp)^{amq^c + arriq^s)] 

rorus [1 {bA - 6p)amq,s] [. 


( 6 . 11 ) 


Here, the flavor indices have been replaced by the particnlar quark flavor (c for 
charm and s for strange). The relation can also be applied to the mass degenerate 
case, where one obtains 


^ [1 + {b^ - bp)amq^^] . 


( 6 , 12 ) 
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On the other hand, it can be extracted directly from the corresponding bare 
subtracted quark mass via the relation 

roM^lrn^^^ = ZMZroiriq^c [1 + &ma"iq,c] • (6.13) 

The continuum extrapolations of all the three dehnitions of the renormaliza¬ 
tion group invariant charm quark mass are illustrated in figure 6.7. Within 
the errors, they all converge to the same continuum result. The hnal result was 
taken from the continuum limit of since the smallest cutoff effects were 

observed in this case. The corresponding linear fit has y^/d.o.f. = 0.32 and the 
extrapolated value at a = 0 is 


roM, = 4.05(7) . 

Using To = 0.5 fm, this translates into =1.597(28) GeV. By inverting equation 
(2.31) with MAPLE, this value can be translated to the MS-scheme and one gets 
^ci^c) = 1.27(3) GeV. In this step, the 4-loop anomalous dimensions and 
f3 have been employed and the error of Aqcd = 238(19) MeV [?] has been taken 
into account. 

The Particle Data group [?] gives md^c) ~ 1.15 — 1.35 GeV as a range for 
the charm quark mass, excluding current data from the lattice. In addition, they 
also quote mc(mc) = 1.26(13)(20)GeV as the world average from the lattice. 
Here, the second error is an estimated uncertainty of 15% due to the quenched 
approximation. 

The determination of the charm quark mass in lattice QCD has a long history 
[?,?,?,?,?,?,?]. Other recent measurements in quenched lattice QGD include 
[?] with Wicijnd) = 1.26(3)(12)GeV and [?] with Wicijnd) = 1.33(8)GeV, with 
statistical and systematic error in the first case and combined statistical and 
systematic error in the second case. The first result has been obtained in the 
0(a)-improved theory at the rather coarse lattice spacing a ~ 0.07 fm and the 
latter has been obtained from a continuum extrapolation. In [?], the continuum 
limit for the charm quark mass was taken by using a step scaling method. The 
final result quoted there is fricijnc) = 1.319(28) GeV. A result for the charm quark 
mass which was obtained in the same way as here, but without a simulation at 
/ds has been carried out in [?], with a final value of Wicijnd} = 1.301(34) GeV. 

An unquenched computation has been carried out by UKQGD [?] with Nf = 3 
flavors of dynamical light quarks. Although only preliminary and at very large 
lattice spacing (a ~ 1.1 fm), no effects of sea quarks compared to the quenched 
results [?,?] were observed. 

The value for the charm quark mass obtained here is compatible within errors 
with all previous determinations. 
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(a/ro)^ 

Figure 6.7: Continuum extrapolation for the renormalization group invariant 
charm quark mass M^. Only the data at the hlled squares entered the respective 
£t. 

6.2.4 Quenched scale ambiguity 

In order to assess the quenched scale ambiguity (cf. section 5.1.1), the data 
has been analyzed again, with the line of constant physics for the continuum 
extrapolation dehned at the physical meson masses = 5.486, to which the 

data has been interpolated. The results and the ambiguity with respect to the 
case tq = 0.5 fm are the following: 


ambiguity 


= 205.8(70) MeV ^ 

-9% 

1^/ 

213(17) MeV ^ 

-11% 

(mo* - 'OT-oJlr' 

= 111.9(66) MeV ^ 

-18% 


= 1.690(33) GeV ^ 

+6% 

mc(mc) 

1.32(3) GeV ^ 

+4% 
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observable 

experiment [?] 

lattice 

precision 

(quenched) 

Fd. 

267(33) MeV 

226(7) MeV 

3% 

Fd* 


239(18) MeV 

8% 

Fds/Fdj 


0.939(61) 

7% 


1.9683(5) GeV 

input 


mo* - mo. 

143.8(4) MeV 

136.0(92) MeV 

7% 

M, 


1.597(28) GeV 

2% 

mc(mc) 


1.27(3) GeV 

2% 


Table 6.6: Summary of results for the Dg- and the D*-meson and the c-quark 
mass. 

The magnitude of the ambiguities agree with equivalent estimates in [?] for the 
decay constant and in [?] for the charm quark mass. It should be mentioned 
here, that the quenched scale ambiguity only gives an estimate for the ambiguity 
inherent to the quenched approximation. The true quenching error can only be 
determined by direct comparison to results in the full theory. 

6.2.5 Discussion 

All results in this section have been summarized in table 6.6 together with the 
hnal error and the corresponding experimental values (where available). Most of 
the observables determined here, are compatible within errors with the previous 
lattice simulations and with QCD sum rules. Especially in the case of the pseudo 
scalar decay constant it has been demonstrated, that it is possible to produce 
lattice data, which matches the precision of experiments, i.e. CLEO-c with a 
predicted error for Fd„ of 2% [?]. 

In the case of the meson mass splitting and the renormalization group in¬ 
variant charm quark mass, the data obtained here is nicely compatible with the 
expected linear scaling in (a/ro)^, which is conhrmed by the small of the 
corresponding hts. 

For the decay constant Fd^, although still acceptable, the relatively large value 
y^/d.o.f. = 1.9 indicates, that the data is not quite compatible with the expected 
linear scaling behavior. In contrast to the mass splitting and the quark mass, the 
boundary-to-boundary correlation functions /J and fy enter the dehnition of the 
pseudo scalar and the vector meson decay constant, respectively. Both quantities 
are known to fluctuate strongly in the Monte-Carlo history. The situation might 
therefore improve with larger statistics at /S^. 

The simulation results from the coarser lattices, which have been excluded 
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from the continuum extrapolation in order to avoid mass dependent cutoff effects, 
deviate considerably from a linear scaling in (a/ro)^. This conhrms the Endings 
in [?] and the upper bound qMq < 0.64 for the heavy quark mass at finite lattice 
spacing, suggested there. 

Despite working in the quenched approximation, most of the results presented 
in this section are compatible within errors with the corresponding values from 
experiment. A slight discrepancy was discovered for the decay constant. However, 
unquenching in this case is expected to increase the value by about 10%. 

Although all sources of systematic error have been taken into account during 
the data analysis, it could be shown, that lattice results can meet the accuracy 
of precision experiments. 


6.3 The Bg- and the B*-meson and HQET 

This section presents the combined analysis of predictions from HQET together 
with the relativistic data for heavy quarks with masses around the charm quark 
mass. 


6.3.1 Incorporation of results from the static limit 

Similar to the observables for the Dg-meson, the quantities 


^3/2 Ips 


R 


Cps’ C'ps/v 


and ro 


Am 


a 


(6.14) 


Spin 


which are defined in section 2.3, have been extrapolated to the continuum at the 
pseudo scalar masses which are collected in table 6.4 together with the associated 
fit range. 

In the case of at the largest mass rupg, the linear continuum extrap¬ 

olation of three data points as depicted in figure 6.8 has a large x^/d.o.f. = 2.5. 
Although P(x^/d.o.f. = 2.5) < 10%, the continuum value has been included into 
the further data analysis, since the extrapolations to a = 0 at all other pseudo 
scalar masses behaved much better. 

For r’Q^^hps/C'ps, the non-perturbatively determined value in the static ap¬ 
proximation [?] 

T'i'ptEGi = 1-74(13) (6.15) 

has been included into the fit as a constraint in the static limit. For the ratio 
R/Cps/y and the mass splitting roAm/Cspm, the asymptotics as detailed in sec¬ 
tion 2.3 have been included to constrain the interpolation. Figure 6.9, 6.10 and 
6.11 show the results. The interpolation has been carried out, once including the 
data in the range of pseudo scalar masses from rompg — r^mp^ and once with the 
masses rgrupg — to^^ps- the larger mass range, a second order polynomial 
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mp 3 with x^/d.o.f. = 2.5. Only the data at the hlled symbols entered the £t. 


was used as the £t-ansatz, while for the smaller mass range also a linear ht was 
applied. The physical point of the Bg-meson and the corresponding results from 
the interpolations are indicated by the circles. Table 6.7 shows the parameters 
of the ht polynomial 

(2i do 

“o “I-1" 7-7? 

romps [rompsy 

The errors associated to the £t-parameter oi in the case of linear interpolations 
is not too large. Thus, Oi can be taken as a rough estimate of the magnitude of 
the hrst-order correction. 

A ht of R/Cps/y over the mass range romp 3 — romp 3 with a 2nd order poly¬ 
nomial is not feasible due to the large error associated with the data. In the 
other cases, 2nd order polynomials ht the data well for both mass ranges. How¬ 
ever, the corresponding coefficients can only be determined very inaccurately and 
significant statements about higher order contributions cannot be made. 


(6.16) 


6.3.2 Interpolation to the Bg -meson 

In the following, the determination of the observables 

Fb,, Fb,/Fbj, andro(mBj - "^bJ (6.17) 

will be presented. One hrst extracts the value of the corresponding interpolation 
at the physical point [?] 

ruB, = 5.3696(24)GeV ^ romp, = 13.604(61), (6.18) 

which is indicated by the blue circles in the respective plots 6.9, 6.10 and 6.11. 
Then, the conversion functions have to be eliminated. In the case of the decay 
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observable 

linear 


quadratic 



Oq tti 

Oq 

ai 

0,2 


mass range mi.. 

. .me 



Cps 

- 

1.7(2) 

-3(2) 

3(4) 

R 

Cps/v 

- 

1.0 

0.4(6) 

-7(3) 

ro Am 

C*spin. 

- 

0.0 

1.8(3) 

1 (1) 


mass range m 3 ., 

. .me 



rl^^Yps 

Cps 

1.7(1) -2.4(7) 

1 .8(2) 

-3(2) 

4(8) 

R 

Cps/v 

1.0 -0.8(2) 

1.0 

- 

- 

ro Am 

C*spin. 

0.0 1.99(8) 

0.0 

1 .6(8) 

2(4) 


Table 6.7: Coefficients and associated error of the polynomial fit oq + romps 
(romps)^ interpolation between the relativistic data and the static approx¬ 

imation for various observables. 
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l/roms, l/^omD, 



l/romB, l/romo. 



1/romps 

Figure 6.9: Interpolation for the decay constant of the pseudo scalar meson with 
a linear and a quadratic £t ansatz (solid and dashed line resp.). Only the data 
at the hlled squares was included into the £t. The circles represent the values at 
the physical point of the Bs-meson. 

















'PS/V ^PS/V 
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1/romB, l/romo. 



l/romo. 
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1/romB, l/romD, 






1/romps 

Figure 6.11: Interpolation for the mass splitting with a linear and a quadratic 
£t ansatz (solid and dashed line resp.). Only the data at the Filed squares was 
included into the £t. The circles represent the values at the physical point of the 
Bs-meson. 
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Fit 

Mass range 

linear 

m^g ... m|s 

quadr. 
m^s ... m|s 

linear 

m^g ... m^g 

quadr. 

m^g ... m|g 


0.50(3) 

0.51(3) 

0.50(3) 

0.53(3) 

Fbs/Fbj 

0.907(9) 

0.99(3) 

0.93(1) 

- 


0.171(5) 

0.16(2) 

0.167(7) 

0.15(4) 


Table 6.8: Results for Bg from the interpolation. 

constant and the ratio of the decay constants, also the square roots of the meson 
masses have to be eliminated (cf. section 2.3). To this end, the conversion 
functions Cx_{MQ/hrg^ were taken at the value of the Yquark mass = 

16.12(29) which has been determined non-perturbatively in quenched QCD [?]. 
The results for all observables for the Bg-meson are summarized in table 6.8 for 
the various interpolations that have been carried out. 

The quadratic interpolation over the mass range mpg... mpg leads to very 
large errors and signihcant results cannot be obtained. The linear interpola¬ 
tion over the larger mass range mpg... m|g on the other hand stretches across 
a domain, where sizeable higher order contributions are to be expected. Both 
interpolations were therefore discarded for the further analysis. In order to de¬ 
termine the hnal results, the linear interpolation including the masses mpg ... mpg 
and the quadratic interpolation including the masses mpg ... mpg were then com¬ 
pared. First, the corresponding results from table 6.8 translated to physical units 
are 


Fit 

Mass range 

linear 

mps • • • 

quadr. 

mpg • • • 

Fb. 

200(10) MeV 

198(9) MeV 

Fbs/Fb* 

0.93(1) 

0.99(3) 


66(3) MeV 

63(6) MeV 


The results of both hts are compatible within la. While the error for the ratio 
of the decay constants and for the mass splitting is twice as large in the case of 
the interpolation with the 2nd order polynomial, the error is approximately the 
same for the decay constant. Comparing the results from both hts, no preference 
becomes apparent. However, in order to arrive at a conservative error estimate, 
the values obtained from the quadratic interpolation were chosen as the hnal 
results. All hnal results for the Bs*^-meson are collected in table 6.10 
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Reference 


FB^^“7MeV Nf 

FB./MeV 

scale 

setting 

ALPHA 

[?] 

206(10) 



ro 

Becirevic et. al. 

[?] 

204(16)(™) 



nip 

Bowler et. al. 

[?] 

220(6)(3) 



U 

CP-PACS 

[?] 

219(10) 

2 

250(10)(1“) 

nip 

CP-PACS 

[?] 

220(4)(31) 

2 

242(9)(!“) 

nip 

de Divitiis 

[?] 

192(6)(4) 



ro 

JLQCD 

[?] 


2 

215(9)(i;S) 

nip 

MILC 

[?] 

199(5)(;i) 

2 

217(6)(;i) 

U 

UKQCD 

[?] 

220(6)(3) 



U 

Wingate et. al. 

[?] 


3 

260(7)(28) 

T 

Ryan [world av.) 

[?] 

200(20) 


230(30) 



Table 6.9: Results for the decay constant Fb, from other groups with statistical 
and systematic errors. 


Using QCD sum rules, Fb„ = 236(30)MeV [?] and Fb^ = 244(21)MeV [?] have 
been determined. Other groups have done lattice simulations in the quenched 
theory as well as with dynamical quarks with Nf = 2 and Nf = 3 flavors. The 
corresponding results have been summarized in table 6.9. Again it turns out, 
that unquenching yields decay constants that are about 10% larger than in the 
quenched case. However, dynamical simulations still suffer from large systematic 
errors. 

The experimental value of the mass splitting is 47.0(26)MeV [?]. 

6.3.3 Quenched scale ambiguity 

As in the case of the Dg- and D*-meson, the quenched scale ambiguity has been 
estimated. The size 


ambiguity 


= 182(10) MeV - 

-8% 

(rrtB* -rrtBjp' 

52(6) MeV - 

-18% 


is roughly the same as for the corresponding observables of the Dg-meson. 
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observable 

experiment [?] 

lattice 

precision 

(quenched) 

Fb. 


198(9) MeV 

5% 

Fbj 


194.664406(1) MeV 

6 % 

Fbs/Fbj 


0.99(3) 

3% 

"iB, 

5.3696(24) GeV 

input 


tub* 

mp* — tuBs 

5.4166(35) GeV 
47.0(26) MeV 

63(6) MeV 

11 % 


Table 6.10: Summary of results for the Bg- and the B*-meson. 

6.3.4 Discussion 

The interpolation between the results from quenched QCD around the charm 
sector and HQET indicates, that the coefficients of the linear and quadratic term 
in the heavy quark expansion 

aoH- — -^7 —+ (6-19) 

romps [rompsY 

are of order 0(1) and one can therefore expect, that HQET is a good approxi¬ 
mation for B(g)-mesons. 

The interpolation has also been used successfully to determine observables of 
the Bg-meson with reasonable errors. All final results have been summarized in 
table 6.10, together with experimental data and the associated errors within the 
quenched approximation. 

The result for the decay constant is compatible with most of the previous 
studies. The mass splitting that was determined from the simulations is not 
compatible with experiment when setting the scale with the Kaon decay constant. 
It agrees however, when setting the scale with the nucleon mass instead. 










Chapter 7 

Summary and outlook 


The work for this thesis was focused on precision measurements of heavy-light 
meson observables in quenched QCD. The systematic errors stemming from 

• discretization effects 

• contributions from excited states 

• the continuum extrapolation 

• finite volume effects 

• interpolation to the physical quark mass 

have been controlled, estimated and considered in the analysis. Only the unknown 
error due to the quenched approximation remains. In particular, the meson decay 
constants and the mass splitting for the Ds - and the Bs -meson and the charm 
quark mass were studied. Moreover, the order of magnitude of the coefficient of 
the leading order contributions to the static approximation in the heavy quark 
expansion was estimated. 

Starting from the MILC-collaboration’s computer program for SU(3) lattice 
gauge theory, a platform independent tool was created and tested, that can ac¬ 
complish all the necessary computations, for example on a PC-cluster, using 
MPI-based parallelism. 

In a scaling study with five lattices of constant volume but decreasing lattice 
spacings a ~ 1, 0.8, 0.7, 0.5, 0.3 fm, the desired observables were extrapolated to 
the continuum. At each lattice spacing, simulations were carried out for six heavy 
quark masses in the region of the charm quark mass, while keeping a seventh 
quark mass at the physical value of the strange quark mass. Deviations from 
the expected linear scaling in (a/ro)^ were observed in all considered observables. 
They were stronger for the heavier meson masses. The expected scaling could be 
recovered by restricting the data that enters the continuum extrapolation with 
the upper bound aMq < 0.64 for the heavy quark mass [?]. 

The results at the physical mass of the Dg-meson show, that a final combined 
statistical and systematic error of 3% on the pseudo scalar meson decay constant 
in the continuum can be achieved. This is at the level of the precision of currently 
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running experiments. For example, CLEO-c claims, that the error on Fd^ from 
experiment will be reduced below 2% in the near future [?]. Also the results for 
the charm quark mass in the continuum could be determined with a hnal error 
of 2%. 

The conversion functions, that are necessary to relate observables in QCD 
and their renormalization group invariant analog in HQET were computed and 
parameterized as functions of the renormalization group invariant quark mass. 
The simulation results of quenched QCD in the region of the charm quark mass, 
extrapolated to the continuum, were successfully combined with predictions from 
HQET by means of an interpolation in the inverse meson mass. Also the error 
due to the finite order in perturbation theory in the conversion functions entered 
the data analysis. 

The results from the interpolation indicate, that the leading spin- and flavor 
symmetry breaking corrections in the heavy quark expansion have coefficients 
that are sufficiently small to expect HQET to be a good approximation for mesons 
containing a 6-quark as the heavy quark. These Endings are compatible with a 
similar study of very heavy relativistic quark masses in small volume [?]. Fur¬ 
thermore, by evaluating the interpolation at the physical point of the Bg-meson, 
predictions for the pseudo scalar and the vector meson decay constant and the 
mass splitting with a combined statistical and systematic error of 5%, 6% and 
10%, respectively, were obtained. 

On a wish list of what has to be done next in order to improve the current 
status are: 

• the calculation of the 1/m-corrections to the static limit 

• the inclusion of a more precise value for the decay constant in the static limit 

• the extension of the calculations to the B-mesons containing a u or d quark 
as the light quark 

• non-perturbative matching 

• the repetition of the simulations with dynamical fermions 

The first two wishes are work in progress by the ALPHA collaboration [?,?] and 
will allow to further constrain the interpolation and thus to reduce the error on 
the observables of the Bg-meson. 

Including the light u- and d-quarks in order to simulate for the B-meson is 
very costly and involves a chiral extrapolation which introduces an additional 
source of systematic errors. 

A program to non-perturbatively match HQET and QCD which would reduce 
the systematic uncertainty due to perturbation theory in the conversion functions 
has been set up by the ALPHA collaboration [?]. 

Finally, only results from simulations of full QCD can be used to reliably 
predict the physical observables, that are necessary for a precision analysis of the 
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Standard Model. When keeping all sources of systematic errors under control, 
this is still a very costly task. 



Appendix A 
Notation 


Pauli matrices 

The Pauli matrices are defined as 

They fulfill the Lie-Algebra 

[cTj, (Jj] ‘2i€ijk(Tk, 


(Ti = 


0 1 
1 0 


, cr2 = 


—t 

0 


, = 


(AT) 


(A.2) 


and are, in context with the iso-spin algebra, also often referred to as r* = Uj. 


Dirac Matrices 

The Dirac Matrices in Euclidean space and in Minkowski space are connected via 


-Euclidean — '-Minkowski 

7i,2,3 — ^ 7i,2,3 


-.Euclidean — -.Minkowski 

7o — 7o 


(A,3) 


They obey the anti-commutation relation 




(A.4) 


and can be constructed from the Pauli matrices. In the chiral representation, the 
Euclidean Dirac matrices 7^ (/r = 0,1, 2, 3) read 


7i,2,3 — 


0 


■*Cri,2,3 

0 


The matrices for /i = 0 and 75 = 7o7i7273 in the chiral representation are 


7o = 


0 1 
1 0 


, 75 = 


1 0 
0 -1 


(A.5) 


(A.6) 


{7m>75} = 


90 


75 anti-commutes with the 7^, 


0 . 


(A.7) 


Appendix B 

The relation between the pole 
mass and the renormalization 
group invariant quark mass 


The pole mass tuq is related to the renormalization group invariant quark mass 
Mq via 


mQ 


mg 

rniWt] 


m[m) 
M, 


Q 


M, 


Q- 


(B.l) 


rn(jn) is the renormalized mass in the MS-scheme of dimensional regularization. 
The ratio mQ/rn{fn) has been determined to three loop precision in [?,?,?] and 
is quoted here for the quenched approximation: 


rriQ 

m(m) 


1 + ^ + 13.4434^ + 190.595|gl 


(B.2) 


The ratio m{m) / Mg on the other hand can be determined from the renormaliza¬ 
tion group equations 


a/i a/i 


(B.3) 


using the 4-loop anomalous dimension of the renormalized coupling [?] and 

the 4-loop quark mass anomalous dimension T^^{g) [?], which have the expansions 


/3MS(^) ^ _ l^^gl _ ^^^9 _ _ _ _ 

T^^{g) = -dog'^ - dig‘^ - d2g^ - d^g^ + ... . 


(B.4) 


The corresponding coefficients are collected in table B.l. Integrating the renor¬ 
malization group equations (B.3), one obtains 


exp 

m[m) 


9(m) 

I ds 


0 


f3^{g) 


do 

bog 



(B.5) 
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[?] 

^MS [7] 

bo 

_ 11 
(47r)^ 

do 

A 00 

to 

hi 

102 

di 

404 

“ (47r)4 

3(47r)4 

b2 

2857 

d2 

2498 

“ 2(47r)6 

(47r)® 

bz 

_ 29243-5033/18 

(47r)® 

dz 

_ 50659 

(47r)^ 


Table B.l: 4-loop anomalous dimension of the coupling and the mass in the 
MS-scheme of dimensional regularization. 


which can be evaluated using MAPLE. 
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Q 

c 

H H 

^ 23 

H ^ 

^ R 

H ^ 
§2 
g N 

hJ 

O < 

gg 

^ s 

m 


b 

X 

HH 

Q 

X 

H 

PLh 

Ph 

c 


non-pert. 

pert. 


parameterization 

parameter error 

reference 


^-loops 



range 


X 


^AiOo) 

_ l-0.8496g^+0.0610gg 

~ l-0.7332gg 


[?] 

X 


Zv{9o) 

_ l-0.7663gg+0.0488gg 

l-0.6369gg 

^<9l<l 

[?] 

X 


ZpW) 

= 0.5233 - 0.0362(/3 - 6.0) + 0.0430(/3 - 6.0)^ 

6.0 < /3 < 6.5 

[?] 

X 


Z{go) 

, n nnncri .1 „2^ l-0.9678ff2+0.04284g4-0.0437336 
- (i + O.OOOOMfiioj l-0.9678gg 

0 < ^0 < 1 0.04% 

[?] 

X 


Zm{P) 

= 1.754(19) + 0.27(10)(/3 - 6) - 0.10(10)(/3 - 6)^ 

6.0 < /3 < 7.0 

cf. (5.5) 


Table C.l: Summary of renormalization constants. 


05 











non-pert. 

pert. 

#-loops 


parameterization 

parameter 

range 

error 

reference 

X 


Csw(fl'o) 

l-0.656c/g-0.152g^-0.054gg 
“ l-0.922gg 

0 < < 1 


[?] 

X 


Ca 

= -0.00756g?nH^ 

0 < < 1 


[?] 

X 


Cv 

- 1 o-oiesssoHSU 

0 < < 1 


[?] 

X 


/9 

6 6.2 6.4 


(stat.)(syst.) 

[?] 


bA 

1.28(3)(4) 1.32(3)(4) 1.31(2)(1) 


X 


bA — bp 

_ -0.00093gg(l+23.3060gg-27.371g;5 

~ l-0.9833gg 

0.881 <gl<l 

< 0.3% 

[?] 

X 


bm 

= (-0.5 - 0.09623^g) 

0 < < 1 

< 1.3% 

[?] 

X 


by 

_ l-0.6518gg-0.1226g;f 

~ l-0.8467c/g 



[?] 


1 

Ct 

= 1 - |0.01346(1)^2 



[?] 


2 

Ct 

= 1 - 0.08900(5)^2 _ 0.0294(3)^^ 



[?] 


Table C.2: Summary of improvement constants. 
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Ki - K 2 

Kl - K 3 

— K 4 

Ki — ^5 

Ki - Kq 

Ki — K 7 


A 

3.07-4.38 

3.26 - 4.56 

3.26 - 4.75 

3.26 - 4.75 

3.26 - 4.75 

3.26 - 4.75 


P2 

3.55 - 4.81 

3.71 - 4.97 

3.71 - 5.13 

3.71 - 5.13 

3.71 - 5.29 

3.71 - 5.29 

romfg{xo) 

Pz 

3.32 - 4.81 

3.32 - 5.08 

3.45 - 5.21 

3.45 - 5.21 

3.32 - 5.21 

3.32 - 5.21 


Pa 

3.13 - 4.28 

3.22 - 4.66 

3.32 - 4.66 

3.32 - 4.86 

3.32 - 4.95 

3.32 - 5.14 



3.34 - 4.22 

3.22 - 4.59 

3.09 - 4.78 

3.09 - 4.84 

3.03 - 4.97 

2.72 - 3.53 


Pi 

3.45 - 4.19 

3.07-4.56 

3.07-4.75 

3.07 - 4.75 

2.89 - 4.75 

2.89 - 4.94 


P2 

3.55 - 4.66 

3.71 - 4.97 

3.55 - 5.13 

3.39 - 5.13 

3.39 - 5.29 

3.39 - 5.44 

rom^{xo) 

Ps 

3.59 - 5.08 

3.45 - 5.21 

3.45 - 5.21 

3.45 - 5.35 

3.32 - 5.35 

3.32 - 5.49 


Pa 

3.61 - 4.66 

3.51 - 4.66 

3.22 - 4.86 

3.22 - 4.86 

3.22 - 4.86 

3.13 - 5.05 


Pb 

3.59 - 4.28 

3.28 - 4.53 

3.16 - 4.66 

3.03 - 4.72 

2.91 - 4.97 

2.66 - 4.72 


Pi 

3.07-4.75 

3.63 - 4.94 

3.82 - 4.94 

4.01 - 5.12 

4.01 - 5.12 

4.19 - 5.12 


P2 

4.02 - 5.13 

4.34 - 5.13 

4.34 - 5.29 

4.50 - 5.13 

4.66 - 5.13 

4.81 - 5.29 

roFps(xo) 

Ps 

3.59 - 5.08 

4.00 - 5.21 

4.27- 5.21 

4.13 - 5.35 

4.27-5.35 

4.40 - 5.21 


Pa 

3.61 - 4.57 

3.89 - 4.76 

4.09 - 4.76 

4.09 - 4.86 

4.38 - 4.95 

4.28 - 5.14 


P^ 

3.34 - 4.53 

3.59 - 4.66 

3.66 - 4.72 

3.66 - 4.78 

3.78 - 4.78 

3.47- 4.41 


Pi 

3.07-4.56 

3.07 - 4.94 

3.07-4.94 

2.89 - 5.12 

2.89 - 5.12 

2.70 - 5.12 


P2 

4.02 - 4.97 

3.39 - 5.13 

3.55 - 5.29 

3.71 - 5.29 

3.87-5.29 

3.87- 5.44 

roFv(xo) 

Ps 

3.18 - 5.21 

4.00 - 5.35 

4.54 - 5.35 

4.54 - 5.49 

4.54 - 5.35 

4.54 - 5.49 


Pa 

3.99 - 4.95 

3.99 - 4.95 

4.18 - 4.76 

4.47 - 4.76 

4.47 - 4.76 

3.80 - 5.05 


P^ 

3.34 - 4.03 

3.28 - 4.22 

3.28 - 4.28 

3.34 - 4.28 

3.28 - 4.47 

3.28 - 4.53 


Table D.l: The plateau ranges in units of tq. 
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Ki — K2 

Ki - K3 

Ki — Ki 

Ki - K5 

Ki - Kq 

Ki — Kj 

a. 

A 

1.556(64) 

1.408(85) 

1.355(98) 

1.32(11) 

1.30(13) 

1.29(15) 

S 

<! 


1.529(56) 

1.396(81) 

1.349(96) 

1.32(12) 

1.29(13) 

1.28(13) 

o 

Pz 

1.571(56) 

1.449(71) 

1.397(87) 

1.36(10) 

1.35(11) 

1.32(13) 


Pa 

1.732(84) 

1.549(99) 

1.48(12) 

1.42(13) 

1.38(15) 

1.37(17) 


Pb 

1.552(92) 

1.50(11) 

1.50(18) 

1.44(17) 

1.37(19) 

1.36(56) 

> 

Pi 

1.35(13) 

1.45(16) 

1.32(18) 

1.37(32) 

1.43(31) 

1.38(25) 

s 

<1 

P2 

1.44(10) 

1.275(72) 

1.30(11) 

1.29(12) 

1.27(14) 

2.667(99) 

o 

P3 

1.431(86) 

1.363(88) 

1.32(11) 

1.30(12) 

1.30(13) 

1.28(15) 


Pa 

1.378(86) 

1.338(99) 

1.41(12) 

1.38(14) 

1.37(16) 

1.34(16) 


P3 

1.338(82) 

1.317(77) 

1.302(94) 

1.29(11) 

1.29(15) 

1.26(24) 

m 

Cb 

Pi 

1.598(34) 

1.424(38) 

1.354(43) 

1.310(47) 

1.282(51) 

1.260(58) 

<] 

P2 

1.534(26) 

1.496(37) 

1.471(40) 

1.428(50) 

1.417(60) 

1.338(62) 

o 

Ps 

1.623(33) 

1.474(32) 

1.420(38) 

1.447(43) 

1.435(49) 

1.390(56) 


Pa 

1.816(30) 

1.692(44) 

1.619(49) 

1.518(52) 

1.521(72) 

1.590(77) 


P3 

1.723(42) 

1.830(56) 

1.804(47) 

1.829(52) 

1.848(69) 

2.233(88) 

> 

fa 

<1 

Pi 

P2 

1.341(61) 

1.89(14) 

1.73(14) 

1.64(15) 

1.61(17) 

1.55(18) 

o 

P3 

2.01(19) 

1.46(12) 

1.240(77) 

1.192(81) 

1.114(73) 

0.973(63) 


Pa 

1.53(11) 

1.514(96) 

1.43(10) 

1.33(15) 

1.35(19) 

1.62(14) 


P3 

1.86(15) 

1.85(16) 

1.89(14) 

1.83(14) 

1.87(14) 

1.91(12) 


Table D.2: Estimates for the gap energy as obtained from fits to the effective 
masses and the decay constants. For roApv at /5i, the data was too noisy to 
allow for a sensible £t. 
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Ki — fi:2 

Ki - K.3 

Kl — 

Kl - Ks 

Ki - Kq 

Ki — K,7 

cn 

A 

3.147(15) 

4.294(19) 

4.971(22) 

5.541(24) 

6.005(26) 

6.439(28) 

pH 

A 

3.069(14) 

4.273(20) 

4.988(23) 

5.597(25) 

6.099(28) 

6.574(30) 

o 

A 

3.093(15) 

4.286(20) 

5.012(23) 

5.636(26) 

6.153(28) 

6.648(30) 


A 

3.095(18) 

4.277(25) 

5.020(29) 

5.678(32) 

6.229(35) 

6.756(38) 


A 

3.742(28) 

5.179(36) 

6.250(42) 

7.170(48) 

8.818(58) 

11.588(77) 

> 

A 

3.645(22) 

4.637(26) 

5.253(29) 

5.766(33) 

6.203(35) 

6.609(38) 

A 

3.567(20) 

4.612(23) 

5.272(26) 

5.843(29) 

6.318(31) 

6.769(35) 

o 

A 

3.604(24) 

4.653(26) 

5.323(29) 

5.908(31) 

6.402(33) 

6.884(36) 


A 

3.637(31) 

4.656(34) 

5.337(37) 

5.954(40) 

6.475(43) 

6.967(45) 


A 

4.141(35) 

5.461(41) 

6.480(46) 

7.370(51) 

8.981(61) 

11.698(77) 

CO 

A 

0.493(12) 

0.528(16) 

0.539(19) 

0.547(22) 

0.551(24) 

0.554(27) 

Ph 

A 

0.524(10) 

0.565(13) 

0.577(16) 

0.583(18) 

0.589(20) 

0.599(23) 

o 

A 

0.539(11) 

0.584(15) 

0.599(18) 

0.607(21) 

0.606(24) 

0.604(26) 


A 

0.5508(98) 

0.594(13) 

0.606(15) 

0.616(17) 

0.628(19) 

0.646(22) 


A 

0.555(11) 

0.572(15) 

0.571(18) 

0.567(20) 

0.558(26) 

0.233(11) 


A 

0.594(36) 

0.565(36) 

0.535(38) 

0.498(42) 

0.471(44) 

0.441(47) 

> 

A 

0.572(26) 

0.578(26) 

0.570(28) 

0.559(31) 

0.551(33) 

0.547(37) 

O 

A 

0.629(44) 

0.663(44) 

0.667(50) 

0.661(53) 

0.648(56) 

0.651(59) 


A 

0.732(49) 

0.700(42) 

0.689(42) 

0.687(45) 

0.687(48) 

0.681(49) 


A 

0.620(39) 

0.600(38) 

0.581(37) 

0.568(39) 

0.550(42) 

0.248(19) 

* w 

Q 

A 

0.829(52) 

0.935(61) 

1.009(73) 

1.097(92) 

1.17(11) 

1.26(13) 


A 

0.917(44) 

0.978(45) 

1.012(49) 

1.043(55) 

1.068(61) 

1.095(68) 

m 

P 

A 

0.857(58) 

0.880(55) 

0.897(61) 

0.919(65) 

0.936(71) 

0.928(71) 

A 

0.753(50) 

0.848(47) 

0.880(48) 

0.897(51) 

0.915(53) 

0.949(55) 


A 

0.895(52) 

0.953(51) 

0.982(53) 

0.999(56) 

1.015(63) 

0.939(63) 

u 

03 

g 

A 

1.657(33) 

3.206(68) 

4.332(96) 

5.43(12) 

6.43(15) 

7.48(18) 


A 

1.578(30) 

3.156(64) 

4.280(89) 

5.35(11) 

6.33(14) 

7.32(16) 

O 

A 

1.657(28) 

3.200(57) 

4.293(79) 

5.33(10) 

6.26(12) 

7.21(14) 


A 

1.691(25) 

3.161(48) 

4.204(65) 

5.193(83) 

6.072(99) 

6.96(12) 


A 

2.394(36) 

4.278(65) 

5.820(91) 

7.23(11) 

9.93(16) 

15.04(26) 


Table D.3: Plateau averaged data for the effective masses, the decay constants 
and the renormalization group invariant quark mass. 
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APPENDIX D. SIMULATION RESULTS 




K,i — K2 

- K 3 

Kl — K4 

Hi — H5 

Hi - Hq 

Hi — H7 

w| CO /9i 

,0hL Cb Ml 
$ (^2 

0.940(21) 

1.053(30) 

1.114(38) 

1.163(45) 

1.200(52) 

1.232(59) 

0.996(19) 

1.127(26) 

1.196(32) 

1.249(38) 

1.295(44) 

1.349(51) 


1 .020(20) 

1.165(29) 

1.243(37) 

1.306(45) 

1.341(52) 

1.370(59) 


Pi 

1.038(19) 

1.185(26) 

1.264(31) 

1.334(37) 

1.402(42) 

1.481(50) 


Pb 

1.082(24) 

1.208(33) 

1.280(41) 

1.335(49) 

1.419(67) 

0.660(32) 


Pi 

0.641(41) 

0.804(54) 

0.894(66) 

0.991(84) 

1.07(10) 

1.16(12) 


1- P2 

0.702(35) 

0.840(39) 

0.896(44) 

0.939(50) 

0.973(56) 

1.005(63) 


^ Ps 

0.661(46) 

0.755(48) 

0.793(55) 

0.826(59) 

0.850(65) 

0.849(66) 


Pi 

0.579(40) 

0.726(42) 

0.776(43) 

0.805(46) 

0.830(49) 

0.869(51) 


P 5 

0.741(44) 

0.845(46) 

0.891(48) 

0.918(52) 

0.946(59) 

0.887(59) 


a Pi 

0.395(15) 

0.333(15) 

0.264(16) 

0.214(19) 

0.181(21) 

0.183(20) 

<1 

cI/92 

0.386(10) 

0.3287(99) 

0.2713(95) 

0.238(10) 

0 .210(12) 

0.193(13) 


0.410(15) 

0.356(15) 

0.300(16) 

0.266(15) 

0.242(16) 

0.233(17) 


Pi 

0.428(19) 

0.367(18) 

0.306(17) 

0.272(18) 

0.242(18) 

0.219(18) 


Pb 

0.389(17) 

0.336(16) 

0.282(16) 

0.252(15) 

0.227(15) 

0.215(15) 


Table D.4: Plateau averaged data for the decay constant, the ratio of the pseudo 
scalar to the vector meson decay constant and the mass splitting converted to 
HQET. 









Appendix E 
Running the code 


This appendix explains the structure of the PC-code, how it has to be compiled 

and how the input files that specify run parameters have to be designed. For the 

official MILC-Code documentation, please refer to 

http : //www. physics. Utah. edu/~detar/niilc/milcv6. html 

or contact the author of this thesis in case of any questions. 

E.l Directory structure 

After unpacking the code, the following directories will be created: 


f_A 


This is the project’s main directory. It contains spe¬ 
cific program code for the computations of correlation 
functions needed in this work. 


schroed_pg 


Gauge-update routines specific for Schrodinger func¬ 
tional boundary conditions. 


generic 


Generic routines, like e.g. I/O or the layout for paraf 
lelization. 


generic_clover Inversion routine (BiGGstab). 


generic_pg 


Generic gauge update routines. 


generic_schroed 


Routines that are specific to Schrodinger Functional 
boundary conditions. 


generic_wilson Routines for boundary sources and the Dirac operator. 
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include 


Macros and declarations and definitions for structures 
like e.g. the site-structure. 


library 


Linear algebra routines. 


E.2 Job steps in a production run 

A production run is organised as follows: There is always one binary for the gauge- 
updates. It starts from a cold gauge conhguration for generating a thermalized 
held conhguration or reloads the conhguration from the previous Monte-Carlo 
step. The resulting held conhguration will be written to disc (In contrast to 
the experience with the APEMille computer at NIC/DESY Zeuthen, the I/O 
is not very time-consuming, even for very large lattices). The binary for the 
calculation of the correlation functions reads the gauge conhgurations from disc, 
computes the propagators and evaluates the correlation functions. If desired, the 
computed propagators can be saved to disk. All run-parameters are handed over 
to the programs by input hies via STDIN. 

E.3 Compiling the code 

Fist a list of important compiler hags: 

DOUBLE: 

If it is dehned in the hie include/config.h, all arithmetics will be done 
with double precision arithmetics. If it is not dehned, only some global sums 
will be done in double precision. Note however, that setting this hag will 
double the CPU-time and will also nearly double the amount of memory to 
be allocated. 

FIELD.MAJOR and TMP.LINKS: 

If dehned in Make_template, the code will be compiled for use of held major 
(cf. section 4.2.3). This improves performance. Additional memory for the 
allocation of temporary helds is needed. 

FORWARD and BACKWARD: 

Depending on which hag is dehned, only the propagators in the forward 
or in the backward direction will be computed. This is convenient for the 
simulation of large lattices with long run-times. At the HLRN, this allowed 
to submit the job-steps into job-queues with a shorter waiting time because 
of reduced wall-clock time. For smaller lattices just set both, the FORWARD 
and BACKWARD hag. Note, that saving propagators works only for one of 
the hags being dehned. 
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Depending on the local installation of the MPICH library, one may have to 
change the PATH to the correpsonding libraries in the Makte_linux_mpi file. 

Comilation of the Pure Gauge part: 


# cd schroed_pg 

# make -f Make_linux_mpi su3_schr_ora 
Compilation of the Inversion routine 

# cd f_A 

# make -f Make_linux_mpi su3_schr_cl_bi 

The directories schroed_pg and f_A should now contain the executables 
su3_schr_ora and su3_schr_cl_bi. 

E.4 Job Scripts 

Sample job script (DESY-Zeuthen-Cluster): 


################################################################## 
#!/bin/csh 

# QSUB -e multiple.err 

# QSUB -r jobtest 

# jump into working-directory 
cd <your working dir.> 

# start the thermalization procedure 
mpirun -np <no. of cpus to use> -machinefile 

<file containing machine names> 

./bin/su3_schr_ora_single \< ./<thermalize input file> » 
thermalize.out 

# Keep thermalized configuration and copy it to working 

# configuration 

cp thermalized checkpoint 

# define a variable that counts job steps 
set a=(l) 

# start loop over job-steps (e.g. 100 measurements) 
while ($a <= 100) 
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# start job 

# Take care for seed, increase it before every update 
awk ’{if($l ~ /iseed/)print $1,$2+123;else print $0}' 

<update input file> > dummyfile 
cp dummyfile <update input file> 

mpirun -np <no. of cpus to use> -machinefile 
<file containing machine names> 

./bin/su3_schr_ora \< ./<update input file> »update.out 

mpirun -np <no. of cpus to use> -machinefile 
<file containing machine names> 

./bin/su3_schr_cl_bi \< <measurement input file> 
»measure. out 


# increase a by one 

set a=(‘expr $a + 1‘) 

end 

rm dummyfile 

################################################################## 

The file <thermalize input file> may be the following (It is important to re¬ 
move all comments!): 


################################################################## 
prompt 0 

nx 20 # define lattice dimensions 

ny 20 
nz 20 
nt 20 

iseed 12318352 # seed for random number gen. 


warms 0 
trajecs 500 
traj_between_meas 5 
beta 7.8439 
bc_flag 0 

steps_per_trajectory 10 

qhb_steps 1 

fresh 

save_serial 
./thermalized 


# warm ups 

# no. of trajectories 

# output gauge-info every # steps 

# beta 

# boundary condition flag 

# heatbath steps 

# OR steps 

# start with flat gauge config 

# save gauge config as binary after end 

# name of gauge-config file 


################################################################## 
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The file <update input file> may be the following (It is important to remove 
all comments!): 


################################################################## 
prompt 0 

nx 20 # define lattice dimensions 

ny 20 
nz 20 
nt 20 

iseed 12318352 # seed for random number gen. 


warms 0 
trajecs 25 

traj_between_meas 25 
beta 7.8439 
bc_flag 0 

steps_per_traj ectory 
qhb_steps 1 
reload_serial 
./checkpoint 
save_serial 
./checkpoint 


# warm ups 

# no. of trajectories 

# output gauge-info every # steps 

# beta 

# boundary condition flag 
10 # heatbath steps 

# OR steps 

# start with flat gauge config 

# save gauge config as binary after end 

# name of gauge-config file 


################################################################## 


The hie <measurement input file> may be the following! (It is important to 
remove all comments!): 


################################################################## 


prompt 0 
nx 20 
ny 20 
nz 20 
nt 20 


# define lattice dimensions 


number_of_kappas 5 # 

bc_flag 0 # 

num_smear 0 # 

# 

kappa 0.133373 # 

kappa 0.128989 
kappa 0.128214 
kappa 0.127656 


total number of hopping params. 
boundary condition flag 
this is redundant, 

I will remove it soon 
values of the hopping params. 
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kappa 0.125309 


cttilde 0.9862 
clov_c 1.3066 
ferni_phases [0] 0.5 
ferm.phases[1] 0.5 
ferm.phases[2] 0.5 
max_cg_iterations 100000 
max_cg_restarts 2 
error_for_propagator le-14 
error_for_propagator le-14 
error_for_propagator le-14 
error_for_propagator le-14 
error_for_propagator le-14 
reload_serial 
<filename> 
num_prop_load aa 
which_prop_load bb 
# numers of the hopping params. 


# cttilde 

# csw 

# theta 


# solver residuals 


# reload binary field config. 

# name of binary field config. 

# number of propagators to reload 

# aa consecutive lines with 


# of the propagator which you 

# want to save 

reload_serial_prop <filename> # conescutive lines with the names 

# of the propagators 

num_prop_sav cc # number of propagators to save 

which_prop_sav bb # same as for load 

save_serial_prop <filename> # 

################################################################## 


If one of the lines num_prop_load or num_prop_sav has 0 inpnt value, leave 
out the which_prop_load/sav and reload/save_serial_prop lines 


E.5 Hints 


• For running the program on a single CPU, one can either use the mpi-code 
and mpirun -np 1 or compile the code as ” vanilla” version with Make_vanilla. 

• Two test input hies for the inverter are located in the directory f_A/ 
(in.cswl.7 und in.test) The output of the corresponding runs can be 
found in the hies data.cswl.7 and data.test. 


• There are some PERL scripts which the user might be interested in; 
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— singdoub: converts a single precision arithmetics gauge config into a 
double precision. 

— alpha2milc: a tool that converts binary or ASCII ALPHA-collaboration 
(APEMille for binary) gauge conhgurations into a format readable by 
the MILC code 

— mouta: extracts correlation function data from output hies and pro¬ 
duces output that is directly readable by the data analysis program 
used here 
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Erratum 


This addendum summarizes the changes and corrections to the thesis (apart from 
typos and corrected references), that were carried out as suggested by the referees. 

Page Comment 

p. 1 The photons of course couple to both, left- and right-handed fermions, 

p. 2 The branching ratio (1.3) has been replaced by its tree-level expres¬ 

sion. 

p. 14 In formula (2.23), the r.h.s. has been corrected. The factor 

{2bQg'^{mQ))'^o' wrong term in the exponent have been 

removed. 

p.l5 There was a wrong sign in formula (2.30) 

p. 17 The definition of the vector meson decay constant (2.38) has been 
corrected. 

p. 28ff The iso-structure of the currents has been formulated more consis¬ 
tently. 

p. 47 The integrated auto-correlation time Tint is now given normalized with 
the total number of intermediate OR and update steps, 
p. 52ff The arguments of the functions Zp have been corrected and rewritten 

in a less misleading way. The error estimates for Zp and Zm have 

been corrected in (5.4), (5.5) and figure (5.1). 
p. 73, The way, the scale has been set in the various lattice simulations for 
p. 85 Fd„ and Fb^ has been added to the tables (6.5) and (6.9). 

p. 80 The values of the errors in table 6.7 have been corrected. 

In addition, reference [?] was added to table 6.9 on page 85. 
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